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SECTION  I 


INTRODUCTION  AND  SUMMARY 


This  report  accompanies  a complete  computer  code  to  directly  calculate 
the  current  density  and  charge  density  induced  on  the  model  of  an  aircraft 
depicted  in  figure  1.  A description  of  the  field  incident  on  the  model 
is  presented  in  figure  2.  The  intent  is  to  calculate  these  induced  den- 
sities for  a continuous  spectrum  that  includes  the  low  and  resonant  range 
of  frequencies.  The  background  establishing  the  need  for  the  direct  cal- 
culation of  these  quantities  for  EMP  external  coupling  purposes  is  presented 
in  reference  1.  In  that  report,  it  was  demonstrated  that  serious  flaws 
existed  in  the  prevalent  stick  model  approach  used  to  perform  EMP  external 
interaction  calculations  for  aircraft.  The  essence  of  the  error  in  that 
approach  was  the  assumption  that  the  current  density  distribution  on  the 
aircraft  could  be  simply  obtained  from  a knowledge  of  the  bulk  current. 

For  this  reason,  the  stick  model  approach  was  intended  only  to  calculate 
the  total  current.  In  reference  1,  it  was  demonstrated  that  methods 
employed  to  obtain  the  current  density  distribution  from  the  bulk  current 
made  their  largest  error  at  low  frequencies  for  which  the  density  calculation 
was  generally  assumed  to  be  the  most  accurate.  The  fact  that  an  error  becomes 
larger  when  it  was  expected  to  become  smaller  is  an  indication  that  the 
physics  of  the  problem  was  not  totally  understood.  Specifically,  it  was 
shown  that  the  stick  model  approach  could  not  yield  a low  frequency  cur- 
rent density  that  was  compatible  with  a magnetostatic  solution. 

The  significance  of  the  fact  that  current  densities  were  incorrectly 
calculated  will  now  be  explained.  It  is  the  current  density  rather  than 
the  bulk  current  that  is  required  by  commonly  employed  deliberate  antenna 


1.  Sancer,  M.  I.,  R.  W.  Latham  and  A.  D.  Varvatsis,  Relationship  Between 
Total  Currents  and  Surface  Current  Densities  Induced  on  Aircraft  and 
Cylinders.  Interaction  Note  194,  Air  Force  Weapons  Laboratory,  August 


1974. 


and  aperture  analyses  for  the  subsequent  calculation  of  voltages  and 
currents  driving  mission  critical  subsystems  contained  within  the  air- 
craft. We  note  that  charge  density,  another  important  external  coupling 
quantity,  was  calculated  by  employing  the  continuity  equation  without 
considering  the  derivative  of  the  azimuthal  component  of  current  density. 

Before  dismissing  the  concept  of  obtaining  useful  external  interaction 
results  from  approaches  based  on  inferring  information  from  a calculated 
bulk  current,  it  is  necessary  to  consider  a new  approach  to  the  problem 
(ref.  2).  In  reference  2,  the  current  density  and  charge  density  are 
calculated  in  the  old  manner  based  on  the  bulk  current;  however,  a 
"correction  solution"  is  added  to  the  old  current  density.  The  addition 
of  this  correction  solution  was  demonstrated  to  improve  the  agreement 
between  the  calculated  current  density  and  that  measured  on  a laboratory 
model.  It  is  not  clear  from  reference  2 what  the  correction  solution  means. 
It  is  essential  that  the  meaning  be  totally  understood  in  order  that  the 
potential  and  the  limitations  of  the  approach  be  understood.  The  correction 
solution  employed  in  reference  2 is  an  approximation  to  the  magnetostatic 
current  density  solution  for  the  aircraft  model.  This  approximation 
corresponds  to  the  magnetostatic  solution  for  an  infinite  circular  cylinder. 
Consequently,  it  can  only  be  expected  to  be  a good  approximation  on  circular 
cylinder  portions  of  the  aircraft  model  that  are  not  in  the  proximity  of 
any  other  aircraft  feature.  One  can  obtain  an  appreciation  for  the  limita- 
tion of  this  approximation  by  noting  that  it  is  not  capable  of  calculating 
the  azimuthal  component  of  the  current  density.  Depending  on  such 
quantities  as  the  incident  field,  frequency,  and  location  of  the  body, 
the  azimuthal  component  of  the  current  density  can  be  larger  than  the 
longitudinal  component. 

By  understanding  that  the  correction  solution  is  an  approximation  to 
the  magnetostatic  solution,  one  can  immediately  appreciate  one  of  the 
hidden  difficulties  of  the  approach.  It  is  as  difficult  to  obtain  a 


2.  Taylor,  C.  D.,  K.  T.  Chen  and  T.  T.  Crow,  Electromagnetic  Pulse 
Interaction  with  the  EC-135  Aircraft,  AFWL-TR-75-205,  Air  Force 
Weapons  Laboratory,  June  1976. 
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magnetostatic  solution  for  an  aircraft  model  that  includes  important 
features, as  it  is  to  use  our  approach  to  calculate  the  current  density 
directly.  In  fact,  we  obtain  the  magnetostatic  current  density  for  our 
aircraft  model  simply  by  setting  the  frequency  equal  to  zero  in  our  code. 
We  see  a possible  benefit  of  adding  our  magnetostatic  current  density  as 
the  correction  solution  in  order  to  obtain  a better  approximation  to  the 
current  density  using  the  modified  bulk  current  approach  of  reference  2. 
The  benefit  of  using  our  result  in  this  manner  would  only  be  to  save 
computer  running  time  since  we  have  to  do  as  much  analysis  and  programming 
to  obtain  the  magnetostatic  solution  as  to  obtain  the  dynamic  solution. 

The  errors  obtained  would  only  be  quantifiable  by  running  our  dynamic 
code  at  selected  frequencies  and  comparing  these  results  to  those  obtained 
by  the  improved  approximate  method.  Finally,  we  note  that  the  addition 
of  the  magnetostatic  solution  cannot  improve  the  calculation  of  the  charge 
density.  This  is  the  case  because  the  divergence  of  the  magnetostatic 
solution  is  zero. 

The  magnetostatic  solution  has  considerably  more  significance  than 
has  already  been  discussed  and  this  is  a topic  dealt  with  in  a recent 
report  (ref.  3).  In  that  report  it  was  demonstrated  that  the  current 
density  induced  on  a metallic  body  by  an  incident  monochromatic  plane  wave 
behaved  predominantly  like  a magnetostatic  solution.  It  was  shown  that 
this  was  the  case  for  an  extended  band  of  low  frequencies  and  this 
frequency  band  can  in  turn  be  shown  to  correspond  to  a significant  portion 
of  the  energy  contained  in  a typical  EMP  spectrum.  It  is  concluded  that 
the  magnetostatic  response  of  a metallic  system  (aircraft,  missile,  etc.) 


3.  Sancer,  M.  I.,  Fundamental  Errors  Associated  with  the  Gross  Modeling 
of  the  Physical  Features  of  Metallic  Enclosures,  AFWL-TR-76-297 , 

Air  Force  Weapons  Laboratory,  December  1976. 
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should  be  considered  In  determining  low  frequency  modeling  requirements 
for  EMP  external  coupling  purposes.  The  primary  significance  of  magneto- 
static related  modeling  requirements  is  that  they  are  much  more  severe 
than  what  has  previously  been  thought  justifiable  as  a result  of  con- 
sidering long  wavelength  far  zone  scattering  results.  An  appreciation 
concerning  the  impact  of  magnetostatics  on  determining  modeling  require- 
ments can  be  obtained  by  noting  the  exact  analogy  between  magnetostatics 
and  irrotatlonal  and  incompressible  fluid  flow  around  a rigid  body.  For 
a rigid  perfectly  conducting  body,  the  velocity  flow  lines  and  the  magnetic 
field  lines  in  the  vicinity  of  the  body  are  identical.  If  certain  features 
of  aircraft  or  missiles  would  have  a significant  effect  on  fluid  flow, 
then  they  would  have  a significant  effect  on  the  current  density  induced 
by  an  EMP.  One  should  become  concerned  with  modeling  such  features  as 
missile  tips  and  fins  as  well  as  aircraft  features  such  as  engines,  wing 
cross  sections,  and  extended  junctions,  particularly  when  the  point  of 
entry  is  in  the  prr  ''ity  of  these  features. 

In  the  process  of  developing  the  computer  code  described  in  this  report, 
we  became  aware  of  the  magnetostatic  modeling  consequences  just  discussed. 
This  fact  had  an  effect  on  our  philosophy  and  approach  in  developing  the 
computer  code. 

First,  wherever  possible  we  developed  general  results  before  restricting 
attention  to  our  particular  model.  Second,  we  realized  that  our  approach 
was  dependent  upon  making  zone  sizes  small  enough  to  approximate  the  surface 
current  density  components  as  constants  over  the  zone.  In  order  to  accom- 
modate the  possibility  of  rapidly  changing  current  densities  in  the  vicinity 
of  edges,  junctions,  and  high  curvature  we  built  a great  deal  of  input 
controllable  zoning  into  the  code.  This  allows  numerical  experimentation 
to  determine  the  effect  of  varying  the  zone  size  in  the  proximity  of  the 
above  features  and  is  potentially  useful  for  obtaining  zoning  information 
for  more  complex  models.  An  additional  benefit  of  the  zoning  flexibility 
is  that  it  permits  one  to  increase  the  density  of  zones  over  the  entire 
model,  at  the  cost  of  increasing  computer  time,  to  either  improve  the 
accuracy  of  the  solution  or  to  obtain  solutions  for  higher  frequencies. 
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A description  of  the  contents  of  the  remaining  part  of  this  report 
will  now  be  presented.  In  Section  II  we  present  the  magnetic  field 
Integral  equation  (MFIE)  for  the  current  density  and  utilize  the 
symmetry  plane  of  the  aircraft  to  transform  this  equation  into  a pair 
of  more  useful  (for  numerical  purposes)  integral  equations  for  suitably 
defined  fictitious  current  densities.  The  numerical  solution  for  these 
integral  equations  only  requires  zoning  of  half  of  the  aircraft.  In 
Section  III  we  trace  the  path  that  led  us  to  the  final  zoning  scheme 
on  the  surface  of  the  aircraft.  Section  IV  compares  numerical  solutions, 
obtained  by  using  our  MFIE  patch  zoning  approach,  with  experimental 
data.  These  data  were  presented  in  two  recent  reports  (refs.  4,5)  and 
even  though  they  refer  to  measurements  of  the  current  density  on  a 
finite  metallic  circular  cylinder,  rather  than  an  aircraft  model,  they 
demonstrate  the  capability  of  our  approach. 

In  Appendix  A we  demonstrate  certain  properties  of  the  surface 
current  density  induced  on  a perfectly  conducting  body  possessing  one 
or  more  symmetry  planes.  These  properties  are  useful  in  that  they 
reduce  computation  time  and  also  give  insight  as  to  the  distribution 
of  the  Induced  current  density  on  a perfectly  conducting  symmetric 
body.  Part  one  of  Appendix  A considers  a body  with  three  planes  of 
symmetry  and  places  no  restriction  on  the  frequency  of  the  incident 
wave.  It  assumes  a wave  vector  perpendicular  to  a plane  of  symmetry 
and  the  electric  field  parallel  to  an  axis  of  symmetry.  The  resulting 
relationships  involve  current  densities  at  points  symmetric  to  planes 
of  symmetry  other  than  the  plane  perpendicular  to  the  wave  vector.  Part 
two  considers  the  magnetostatic  limit  for  a body  with  only  one  symmetry 


4.  Burton,  R.  W.,  R.  W.  P.  King  and  D.  Blejer,  Surface  Currents  and 
Charges  on  a Thick  Conducting  Tube  in  an  E-Polarized  Plane-Wave  Field, 
II.  Measurements,  progress  report  on  contract  F29601-75-C-0019, 
AFWL/ELPE,  Kirtland  Air  Force  Base,  New  Mexico,  1976. 

5.  King,  R.  W.  P.,  Surface  Currents  and  Charges  on  a Thick  Conducting 
Tube  in  an  E-Polarlzed  Plane-Wave  Field,  IV.  Generalization  to 
Cylinders  of  Various  Lengths,  progress  report  on  contract  F29601-75- 
C-0019,  AFWL/ELPE,  Kirtland  Air  Force  Base,  New  Mexico,  1976. 
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plane  and  an  Incident  plane  wave  with  a wave  vector  perpendicular  to 
the  symmetry  plane.  We  show  that  the  surface  current  density  at  points 
symmetric  to  the  plane  of  symmetry  are  simply  related.  Since  magneto- 
statics is  a good  approximation  for  an  appreciable  band  of  low  frequencies, 
it  follows  that  these  relationships  provide  insight  into  the  distribution 
of  the  current  density  on  the  surface  of  structures  with  a plane  of 
symmetry  including  aircraft. 

Appendix  B presents  the  calculation  of  the  surface  current  density 
induced  on  a perfectly  conducting  ellipsoid  immersed  in  a magnetostatic 
field.  This  calculation  was  very  helpful  in  predicting  the  magnetostatic 
current  distribution  on  a finite  elliptical  cylinder  which  in  turn  was 
used  to  test  the  numerical  results  of  the  various  zoning  schemes  as  we 
explain  in  Section  III. 

Appendix  C presents  the  method  for  numerical  solution  in  detail. 
Specific  topics  treated  are  the  model,  zoning  scheme,  matrix  equations 
for  the  current  density,  coordinates  of  centers  and  boundaries  of  zones, 
calculation  of  matrix  elements,  self-zone  interaction  considerations, 
interpolation  scheme,  calculation  of  charge  density  as  well  as  edge 
and  junction  behavior.  Finally,  Appendix  D supplements  Appendix  C 
by  giving  the  coordinates  and  matrix  elements  for  zones  adjacent  to 
junctions. 


SECTION  II 


FORMULATION  OF  THE  PROBLEM 

The  model  of  the  aircraft  we  employ  is  depicted  in  figure  1.  All 
components  shown  (fuselage,  wings,  horizontal  stabilizers,  vertical  stabil- 
izer) are  perfectly  conducting  elliptical  cylinders  making  perfect  electric 
contact  at  the  intersections.  The  aircraft  is  illuminated  by  a monochromatic 
plane  electromagnetic  wave  of  arbitrary  direction  and  polarization.  We 
are  interested  in  calculating  the  surface  current  and  charge  densities 
everywhere  on  the  surface  of  the  aircraft.  To  do  so  we  employ  the  magnetic 
field  integral  equation  for  the  current  density  and  utilize  the  symmetry 
of  the  aircraft  about  the  xz-plane  to  transform  this  equation  into  a pair 
of  more  useful  (for  numerical  purposes)  integral  equations  for  suitably 
defined  fictitious  current  densities.  Each  of  these  equations  lends  itself 
to  a numerical  solution  with  a matrix  (N/2)  x (N/2)  where  (N  x N)  is  the 
matrix  of  the  original  equation. 

The  magnetic  field  integral  equation  for  a perfectly  conducting  body 
is 

J(r)  - J^(r)  + f K(r irj  • J/r^)  dSQ  (1) 

JS 

where  J^(r)  is  the  incident  current  density,  which  serves  as  the  source 
for  the  integral  equation,  given  by 

J^Cr)  - n(r)  x H^r)  (2) 

n(£)  is  the  unit  normal  to  the  surface  at  £,  H^(r)  is  the  incident  magnetic 
field  given  by 


Site)  - • 


ik  *r 


(3) 
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is  the  wave  vector,  K is  the  kernel  given  by 

K(r;io)  - n(r)  x (VG^r;^)  x I]  - n(r)  x [R  x I]  Q(R)  (4) 

Go  " exPUk0R]/^^R  is  the  free  space  Green's  function,  ^ is  the  identity 
operator  and 

ik  R 
e 0 

Q(R)  - r-  (ik  R - 1) 

47TR3  0 

R ■ r - r 
— — — o 

R - | R|  (5) 


The  surface  S extends  over  all  components  of  the  model  shown  in  figure  1 
including  the  flat  caps. 

It  has  been  shown  in  reference  6 that  for  bodies  possessing  a symmetry 
plane  say  the  xz-plane,  equation  (1)  can  be  transformed  into  the  following 
pair  of  equations 


\ J+(r+)  - J^(r+)  + f I+  • 1 


J+(r  ) dS0 


yj  (r  ) - (r+)  + J K-  (_r+  ;r^)  • j" 


(r  ) dSn 


(6) 


6.  Sancer,  M.  I.  and  A.  D.  Varvatsis,  Analytical  and  Numerical  EMP 
Coupling  Solutions  for  A Class  of  Structures  Attached  to  the  Wing 
of  an  Aircraft.  AFWL-TR- 74-298,  July  1975  (also  published  as  AFWL 
Interaction  Note  197,  October  1974). 
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where  J±(r+)  - [J(r+)  ± • J(Ry  • r+)] 


(7) 


^L)  * \ [^(JL)  ± Ry  * iiCSy  • I>1 


(8) 


I»-/+  +\  . , + + . . + +. 

K (r  ;r  ) - A(r  ;r  ) ± B(r  ;r  ) 
o * o — — o 


(9) 


. , + + . + +. 
;i„)  - £(£  ;io> 


+ +v 


Btr'jrT)  - K(r+;R  • r+)  • R 

— — — o — — — y — o —y 

Jy  is  a reflection  operator  such  that 


(10) 


R - I - 2e  e 

=y  - y y 


JRy  • £.+  * (x,-y,z),  r+(x,y^0,z) 


R • a * (a  ,-a  ,a  ),  a ■ (a  ,a  ,a  ) 
”7  — x y z — x y z 


(ID 


and  S+  is  the  surface  of  the  body  that  corresponds  to  y _>  0.  Thus 
equations  (6)  are  defined  over  the  y > 0 half  of  the  aircraft.  Once  we 

t + 

solve  for  jJ  at  £ (x,y ^0,z)  we  can  use  equation  (7)  to  calculate  at 

jr+  and  R • jr+,  i.e.,  everywhere  on  the  surface  of  the  body: 


J(r+)  - J+(r+)  + J"(r+) 


J(r')  - Ry  • (J+(r+)  - J"(_r+)] 


(12) 


where 


L ■ £y  * l “ (x.-y.*) • 


We  solve  equation  (6)  for  two  orthogonal  surface  components  s • J.~(£+) 

A i + a 

and  t • J,  (£  ) where  s and  t are  unit  surface  vectors  forming  an  orthonormal 
triad  with  the  unit  normal  n:  n - s x t.  For  each  component  of  the  aircraft 
the  s vector  is  defined  as  the  tangent  unit  vector  at  the  intersection  of 
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a cylinder  with  a plane  perpendicular  to  its  axis,  i.e., 


3r/3<J> 

3 “ | 3jr / 3 | (13) 

This  definition  is  true  for  the  fuselage,  vertical  stabilizer  and  the  y > 0 
wing  and  horizontal  stabilizer.  For  the  y < 0 wing  and  horizontal  stabilizer 
s is  given  by  equation  (13)  with  a minus  sign.  The  t-vector  is  obtained 

A 

through  the  relationship  t » n x s.  Referring  to  figures  1 and  3 we  can  show  that 
Fuselage 


t - 1 
x 


y ■ b^r  sin<|>  z ■ -a^r  cos$ 


cos4> 

’ 0 sy  ' 


t ■ 0 

y 


a^  sin# 
5z  ” t^C*) 


t2  - 0 


r * 1 on  walls 
0 r £ 1 on  caps 


on  walls  and 
caps 


on  walls 


t * 0 
x 


t - 0 
x 


sin$ 

cy  " '^(o) 


sini|i  bj,  cos< 

Nl(«)  Cz  " ~N 


~b^  cosq> 
Cz  ’ ^(<0) 


on  x =*  0 cap 


on  x =■  cap 


Wing  for  y > 0 


x ■ a2r  cos(J>+x02  y - y z ■ -b2r  sin4> 


r - 1 on  walls 
0 _<  r 1 on  cap 


a2  sin$ 

sx  “ " N_($)  Sy  " 0 


b2  cos<)i 

n2($) 


on  walls  and  cap 


t - 0 
x 


t - 1 

y 


t ■ o 
z 


on  walls 


b j cos<t> 

'« 5^T  cy  " 0 


a2  sin$ 


z N2(«) 


on  cap 


Figure  3:  Definition  of  the  elliptical  angle  $ such  that  x - a cos<J), 

y - b sin$.  Notice  that  equal  increments  in  <\>  over  a 
quadrant  do  not  correspond  to  equal  arc  lengths. 
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Wing  for  y < 0 


a2  sln$ 
sx  * 1 m$7” 


s “0 

y 


b2  cos4> 
Jz  " N2(0) 


on  walls  and  cap 


t - 0 
x 


t - -1 

y 


cz  - 0 


on  walls 


b0  cos<}> 
Cx  “ N2(<(>) 


a2  sin4> 

ty  " ° cz  N„(<J>) 


on  cap 


(15) 


Horizontal  Stabilizer 


Same  as  for  the  wing  if  we  change  a0,  b2>  xQ2  to  a^,  b3>  x^. 


Vertical  Stabilizer 


x ■ a^r  cos<J)+Xq^  y ■ b^r  sim}>  z = z 


r » 1 on  walls 
0 r ^ 1 on  caps 


a^  sim|i  b^  cos$ 

sx  ‘ " n4($)  sy  " n4(«»  sx  " 0 


on  walls  and  cap 


t - 0 
x 


t - 0 

y 


t - l 
z 


on  walls 


>b4  cos<J> 
Cx  “ N4(<j>) 


a^  sin<f> 

cy  “ n4(*) 


t - o 
z 


on  cap 


(16) 


2 2 2 2 1/2 
where  N^($)  ■ (b^  cos  <>  + a^  sin  $) 

N2(W  ■ (&2  + b2  cos^)^^ 


N4(*)  ■ (a4  sin^$  + b4  cos^)1^ 


(17) 


Equations  (12)  can  be  rewritten  in  component  form 


Jt(r+)  - J*(jT)  + J‘(r+) 


+ ,_+v  t”/_+\ 


Jt(r  ) - t(r  ) • ^ • [J  ( r)  - J ( r ) ] 


(18) 
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(19) 


V/>  - •£(!  ) + Jl(r+) 


* s(r")  • Ry  • [ J+ (r+)  - J"(r+) 


From  relationships  (14)  through  (16)  we  see  that 
Fuselage 


sy(jL+>  “ sy(£_)»  sz(l+)  m -s2(r_) 


cx(X  ) “ tx(r") 


ty(E+)  - "ty(r+),  t2(r+)  - tz(r") 
Wings  and  Horizontal  Stabilizer 

s (JL+)  ■ -sxCr“),  s (r+)  - -s  (r~) 


on  walls  and  caps 


on  walls 


on  caps 


(20) 


on  walls  and  caps 


ty(l  ) - -ty(r") 

^(L)  " t (L),  Cfr+)  - t ( zf) 


on  walls 


(21) 


Vertical  Stabilizer 


sx(!+)  “ -sx(r_),  sy(r+)  - Sy(r“) 


on  walls  and  caps 


’ ‘.of) 

on  walls 

ft 

X 

ft 

1 

ft 

", 

on  cap 

(22) 

relationships  we  see  that 

in  general 

8x(i+)  " *9x^‘) 

txte+>  ’ ‘xOf) 

Sy (r+)  - 8y(r") 

ty (I+)  - -ty(r") 

»z(r+)  - -s2(r") 

tz(r+)  - t2(r") 

(23) 
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and  consequently 


t(r")  • Ry  - [t^jf)  + ty(r_  ) + t^r  ) ] • - t^r  ) - t^r  ) 

+ + ty(r+)  + t^(r+ ) - t(r+) 

Similarly 

s(r  ) • R - -s(r+) 

Thus  equations  (18)  and  (19)  can  be  rewritten  as 

Jt(r+)  - J+(r+)  + J~ (r+ ) 

Jt(r_)  > J+(r+)  - J~(r+) 

Js(£+)  “ ^(£+)  + Js(£+) 

J,(r“>  “ -[Jg(£+)  - Jg(r+)]  (24) 

In  order  to  solve  numerical  equations  (6)  we  first  write  them  in 

component  form  and  then  transform  them  into  a system  of  simultaneous 

algebraic  equations  for  the  components  of  the  fictitious  current  densities 

evaluated  at  the  centers  of  zones  into  which  we  have  divided  the  entire 

■¥  +■ 

aircraft.  That  is,  we  assume  that  the  current  components  J',  J are 

c s 

constant  over  a zone  and  equal  to  their  values  at  the  center  of  the  zone. 
Appendix  C presents  the  numerical  solution  in  detail. 
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SECTION  III 


HIGHLIGHTS  OF  THE  NUMERICAL  SOLUTION 


In  this  section  we  present  some  of  the  highlights  of  the  path  that 
led  us  to  our  final  computer  code.  In  particular,  we  discuss  how  we 
decided  on  the  number  and  size  of  zones,  and  the  accuracy  we  employ  for 
the  calculation  of  matrix  elements. 

Our  first  task  was  to  decide  on  the  number  of  zones  on  the  surface  of 
our  model:  how  many  we  needed,  what  their  optimum  size  was  as  a function 

of  location  and  whether  the  resulting  matrix  to  be  inverted  would  be  too 
large  to  assure  confidence  that  the  round  off  error  would  be  negligible. 

To  assess  this  error  we  considered  the  problem  of  calculating  the 

current  densities  on  a circular  cylinder,  illuminated  by  a plane  wave 

with  broadside  incidence  and  the  electric  field  polarized  along  the  axis 

of  the  cylinder,  by  employing  two  debugged  computer  codes.  One  took 

+ 

advantage  of  one  symmetry  plane  and  solved  for  J~  defined  in  Section  I. 

The  other  utilized  three  symmetry  planes.  The  details  of  the  latter 
approach  are  given  in  Appendix  A where  it  is  shown  that  one  need  only 

_ j | 

calculate  two  fictitious  quantities  and  over  1/8  of  the  surface 

of  the  cylinder  whereas  the  first  approach  requires  calculation  of  J+ 
and  J_  over  half  the  surface.  The  three-symmetry-plane  code  involved 
the  inversion  of  a matrix  108  x 108  whereas  the  one-symmetry-plane  code 
required  the  inversion  of  a (4  x 108)  x (4  x 108)  - 432  x 432  matrix. 

We  displayed  our  results  with  eight  significant  figures  and  did  not 
observe  any  roundoff  error. 

Once  we  gained  confidence  that  large  matrices  of  the  type  generated 
by  our  approach  could  be  Inverted  accurately  we  proceeded  to  determine 
the  minimum  number  of  zones  on  the  aircraft.  To  that  end  we  had  to 
decide  how  many  zones  we  should  use  on  the  elliptical  cylinders  modeling 
the  aircraft  components.  The  choice  of  the  size  of  a zone  is  based  on  the 
requirement  that  the  current  density  components  should  not  vary  appreciably 
over  a zone.  The  variation  of  the  current  density  depends  on  geometrical 
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and  wavelength  considerations.  Geometrical  considerations  determine 
the  minimum  number  of  zones  and  wavelength  considerations  can  increase 
this  number.  The  geometric  requirements  are  exactly  the  ones  required 
in  order  to  obtain  a magnetostatic  solution.  As  an  example,  in  the 
past,  it  has  been  common  practice  to  require  that  the  size  of  a zone, 
both  longitudinal  and  azimuthal,  on  the  surface  of  a circular  cylinder 
should  not  be  larger  than  the  diameter  of  the  cross  section,  independently 
of  the  wavelength.  The  wavelength  requirement  was  that  we  should  have 
at  least  ten  zones  per  wavelength.  Thus  the  wavelength  requirement  was 
automatically  satisfied  via  the  geometry  requirement  for  wavelengths 
larger  than  ten  times  the  diameter  of  the  cylinder.  For  smaller  wave- 
lengths the  wavelength  requirement  determined  the  total  number  of  zones. 

It  should  be  noted  that  the  factor  of  ten  associated  with  the  wavelength 
requirement  is  somewhat  arbitrary.  There  is  no  restriction  built  into 
the  code  based  on  the  factor  being  equal  to  ten.  To  the  contrary,  the 
zoning  flexibility  built  into  the  code  allows  one  to  assess  whether  a 
smaller  factor  (corresponding  to  higher  frequencies  for  fewer  zones) 
yields  sufficient  accuracy. 

For  an  elliptical  cylinder  a geometrical  condition  that  has  been 
imposed  in  the  past  requires  that  the  linear  size  of  a zone  should  not 
exceed  the  minor  axis  of  the  ellipse.  In  order  to  test  this  condition 
we  decided  to  compare  our  computer  code  (with  a variable  number  of  zones) 
to  the  exact  magnetostatic  solution  for  an  ellipsoid  immersed  in  a 
uniform  magnetic  field  (see  Appendix  B) . The  basis  for  the  comparison 
was  that  for  ellipsoids  with  large  c/a  (c  > a > b)  the  variation  of  the 
current  density  around  the  central  cross  section  was  insensitive  to  c/a. 
Thus  we  had  reason  to  believe  that  the  variation  of  the  current  density 
at  the  central  cross  section  of  the  ellipsoid  approximated,  to  a high 
degree  of  accuracy,  the  variation  of  the  current  density  at  the  central 
cross  section  of  a finite  elliptical  cylinder  with  h/a  = c/a  (h  is  the 
cylinder  half  length).  The  computer  code  we  used  employed  all  three 
planes  of  symmetry  for  broadside  incidence  with  an  electric  field  parallel 
to  the  axis  of  the  cylinder.  This  was  done  to  minimize  the  cost  for  each 
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run  as  we  varied  the  number  of  zones  and  naturally  we  ran  the  code  for 
0.  By  examining  the  exact  solution  we  noticed  that  the  current 
density  varied  most  rapidly  near  points  of  high  curvature.  This  meant 
that  our  zoning  should  be  nonuniform  with  zone  size  diminishing  as  we 
approached  the  point  of  highest  curvature.  This  was  an  important 
observation  but  it  led  to  the  following  problem.  In  the  past  we  had 
subzoned  each  zone  in  order  to  accurately  calculate  the  neighboring 
zone  interaction.  As  for  the  salf-zone  calculation  we  had  realized 
the  need  to  analytically  treat  the  integrable  singularity  of  the 
Integrand  in  order  to  secure  a high  degree  of  accuracy.  We  had  accomplished 
chat  by  dividing  the  two-dimensional  integral  into  two  integrals: 
one  with  a two-dimensional  nonsingular  integrand  that  was  calculated 
by  subzoning  and  one  with  a singular  but  integrable  integrand  that 
could  be  reduced  to  a one-dimensional  integral  with  a nonsingular 
integrand.  No  such  precaution  had  been  taken  for  neighboring  zone 
interaction  because  the  integrands  involved  did  not  vary  appreciably 
and  allowed  accurate  calculation  through  subzoning.  For  an  elliptical 
cylinder  with  a large  a/b  (say  larger  than  4)  the  required  zone  size 
in  regions  of  high  curvature  may  be  so  small  that  neighboring  zone 
interaction  could  require  a large  number  of  subzones  to  evaluate  the 
integrals  whose  integrands  now  can  vary  appreciably.  The  code  automatic- 
ally prescribes  the  same  number  of  subzones  for  all  zones  and  this  is 
clearly  unnecessary  for  neighboring  zone  interaction  (in  regions  of  small 
curvature)  where  the  zone  size  is  not  small.  This  would  cause  the  code 
to  be  unnecessarily  costly.  We  decided  to  examine  whether  we  could 
calculate  the  matrix  element  without  sub zoning  the  zones.  We  found 
that  the  two-dimensional  integrals  over  a zone  on  the  walls  or  the  caps 
of  the  cylinder  could  be  reduced  to  one-dimensional  Integrals  with 
nonsingular  Integrands  without  dividing  the  original  integrals  into 
two  Integrals  as  we  did  in  the  past.  (See  Appendix  C,  Section  5 for 
details.)  Thus  we  were  in  a position  to  calculate  the  matrix  elements 
for  self-zone  and  zone-to-zone  interactions  by  evaluating  well-behaved 
one-dimensional  integrals  with  a high  degree  of  accuracy.  An  additional 
benefit  of  making  the  matrix  elements  to  this  degree  of  accuracy  was 


25 


also  achieved.  We  found  that  the  improved  accuracy  of  the  evaluation 
of  the  matrix  elements  improved  the  accuracy  of  the  phase  of  the  solu- 
tion. We  observed  situations  where  it  could  have  been  erroneously 
concluded  that  the  matrix  elements  were  calculated  with  sufficient 
accuracy  based  on  comparing  the  magnitude  of  the  resulting  solution 
with  known  results.  We  calculated  the  same  quantities  with  a code  that 
utilized  more  accurate  matrix  elements  and  noticed  a significant  effect 
on  the  phase  of  the  solution  while  the  magnitude  was  minimally  affected. 

Now  we  will  describe  the  flexibility  built  into  our  code  and  present 
guidelines  for  utilizing  this  flexibility.  A general  description  of  the 
flexibility  is  that  the  zoning  of  the  aircraft  model  can  be  specified 
with  input  data  cards  to  cause  more  dense  zoning  in  regions  near  junctions 
and  edges  as  well  as  the  regions  of  rapidly  varying  curvature  near  the 
leading  and  trailing  edges  of  the  elliptic  cylinder  components.  The 
zoning  density  over  the  entire  aircraft  model  can  also  be  increased 
through  input  data  cards  to  accommodate  frequencies  higher  than  those 
determined  by  the  geometry  limited  zoning.  The  primary  reason  *-hat 
we  have  this  flexibility  is  that  our  method  of  accurately  calculating 
the  matrix  elements  is  insensitive  to  the  size  of  our  zones.  A penalty 
that  is  paid  for  having  this  flexibility  of  nonuniform  zoning  is  that 
the  benefit  of  using  symmetries  to  reduce  the  matrix  generation  time  is 
reduced. 

We  conclude  this  section  by  presenting  guidelines  for  the  geometry 
limited  zoning.  The  presentation  of  these  guidelines  is  facilitated  by 
considering  the  following  categories  of  surfaces  that  require  zoning: 
cylinder  walls,  junctions,  edges  and  end  caps. 

The  zoning  of  the  cylinder  walls  is  first  determined  without  regard 
to  junctions  and  edges.  A description  of  the  initial  zoning  of  the  walls 
is  assisted  by  considering  figure  3.  We  will  describe  the  zoning  of  the 
elliptic  cylinder  walls  with  the  zoning  of  the  circular  cylinder  walls 
being  a special  case.  Our  code  allows  that  the  full  range  for  4>  can  be 
subdivided  as  finely  as  desired.  This  has  utility  for  zoning  in  the 
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vicinity  of  junctions;  however,  for  determining  the  number  and  size  of 
the  zones  independently  of  the  junctions, we  can  consider  that  each  wall 
is  divided  into  either  two  or  four  subdivisions  corresponding  to  a $ 
range  of  90°  (quadrants).  On  the  fuselage  and  on  the  vertical  stabilizer 
we  have  two  quadrants;  on  the  wings  as  well  as  on  the  horizontal  stabilizers 
we  have  four  quadrants.  For  each  of  the  quadrants  we  choose  the  number  of 
<(>  zone  divisions  according  to  the  ratio  of  the  major  to  minor  axis  of  the 
ellipse.  Specifically  the  number  of  zone  divisions  is  the  nearest  integer 
to  half  of  the  sum  of  one  plus  this  ratio.  Once  this  number  is  specified 
as  an  input  parameter,  the  code  chooses  the  size  of  these  zones  so  that 
they  are  smallest  in  the  region  of  the  most  rapidly  varying  curvature. 
Specifically,  it  chooses  them  to  be  equal  increments  in  the  elliptic 
angle  <}>.  The  longitudinal  length  of  a zone  is  taken  to  be  no  larger  than 
the  major  axis  of  the  ellipse.  These  guidelines  are  meant  to  represent 
an  initial  estimate  of  the  geometry  limited  zoning.  The  adequacy  of  the 
zoning  can  be  determined  by  running  the  code  at  zero  frequency  with  the 
initial  geometry  limited  zoning  estimate  and  then  running  the  code  with 
a more  dense  zoning. 

The  zoning  of  junctions,  edges  and  end  caps  is  determined  by  utilizing 
the  zoning  flexibility  of  the  code  and  performing  numerical  experiments. 

We  single  out  the  zoning  of  the  end  caps  for  special  consideration  only 
because  they  are  necessarily  in  the  vicinity  of  an  edge.  Now  we  make 
a further  distinction  between  the  numerical  experiments  in  the  vicinity 
of  junctions  and  edges.  As  explained  in  Appendix  C,  Section  9,  neither 
component  of  current  density  becomes  unbounded  in  the  vicinity  of  a 
junction;  however,  the  component  parallel  to  an  edge  does  become  unbounded. 
We  believe  that  this  singular  behavior  should  be  given  special  analytic 
and  numerical  treatment.  We  attempted  a numerical  subtractive  procedure 
to  treat  this  edge  difficulty  and  obtained  unsatisfactory  results.  Due 
to  time  limitations,  we  have  not  yet  numerically  determined  whether 
a multiplicative  procedure  would  allow  us  to  trust  our  solution  for 
points  arbitrarily  close  to  the  edge.  Despite  the  fact  we  have  misgivings 
aoout  our  solution  arbitrarily  close  to  an  edge,  the  comparison  between 
our  calculations  and  measured  data  (presented  in  the  next  section)  is 
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quite  good  at  a distance  corresponding  to  half  of  an  ordinary  zone 
length. 
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SECTION  IV 


COMPARISON  WITH  EXPERIMENTAL  DATA 


In  this  section  we  compare  numerical  solutions,  obtained  by  using 
our  MFIE  patch  zoning  approach,  with  experimental  data.  Even  though 
these  measurements  and  calculations  were  not  for  the  complete  geometry 
of  the  aircraft  model  they  demonstrate  the  capability  of  the  approach. 

It  should  be  noted  that  the  quantities  that  will  be  compared  correspond 
to  measurements  and  calculations  on  the  surface  of  an  object.  It  is  more 
difficult  to  calculate  surface  distributions  than  far  zone  scattering 
results,  thus  the  excellent  degree  of  agreement  to  be  demonstrated  is 
very  encouraging.  To  further  relate  the  calculations  used  in  the  comparison 
to  our  final  code,  we  note  that  through  the  use  of  input  data  cards  our 
complete  aircraft  code  can  be  reduced  to  yield  results  for  the  geometry 
used  in  the  comparison.  At  this  time  an  experimental  program  is  underway 
to  measure  the  current  density  induced  on  the  complete  aircraft  model  for 
which  we  developed  our  code. 

The  comparison  of  the  calculated  and  measured  data  is  facilitated  by 
considering  figures  4 through  8.  All  of  these  figures  contain  material 
that  was  presented  in  references  4 and  5.  The  curves  in  figures  5,  6a, 

7a,  and  8a  were  traced  from  Xerox  copies  of  curves  presented  in  those 
reports.  The  figure  containing  the  description  of  the  experiment 
pertaining  to  the  data  presented  in  these  figures  is  also  redrawn  based 
on  a figure  presented  in  reference  4.  All  of  these  figures  were  originally 
redrawn  for  presentation  in  reference  3. 

The  intent  of  the  experiment  is  to  simulate  a monochromatic  plane 
wave  incident  on  a tube  having  a total  length  of  2h.  By  referring  to 
figure  4,  we  can  see  that  the  angle  6 is  defined  so  that  0°  corresponds 
to  the  deep  shadow  region  and  180*  corresponds  to  direct  illumination. 

The  6 in  figures  6a,  7a,  and  8a  correspond  to  this  definition  and  z 
is  the  axial  distance  ranging  from  0 at  the  ground  plane  to  h at  the  top 
of  the  cylinder.  The  quantities  | Kz | • |Kq|,  and  6z  plotted  in  these 
figures  are  the  magnitude  of  the  axial  component  of  the  current  density. 
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Figure  5:  Measured  magnitude  and  phase  of  surface  density  of  outside 
axial  current  on  tubular  cylinder  with  open  end,  flat  and 
hemispherical  end  caps.  E-polarization. 
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Figure  6a:  Measured  amplitude  of  axial  surface  density  of  outside 

current  on  tubular  cylinder.  E-polarization  (large 
outdoor  ground  screen) . 


Figure  6b:  Calculated  airplitude  of  axial  surface  current  density 

on  a finite  circular  cylinder. 
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ure  7a:  Measured  amplitude  of  surface  density  of  transverse  current 
on  tubular  cylinder.  E-polarization;  normal  incidence 
(large  outdoor  ground  screen). 


i = 34 

,-  / \ - 

/ 30  \ 

/ /'Xn  \ 

i i i i 

0 30  60  90  120  ISO  180 

9 (deg) 


0 • (deg) 

80x/60  - ,.5 
x/ioo 

Hr 

.// b ...  - 1.0 


y//i2° 
jy// 1-»° 

j:  20 

i60 

«•  *■*  / 


12  18  24 

DISTANCE  Z (cm) 


Figure  7b:  Calculated  amplitude  of  transverse  surface  current  density 
on  a finite  circular  cylinder. 
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the  magnitude  of  the  transverse  component  of  the  current  density,  and  the 
phase  of  the  axial  component  of  the  current  density.  The  same  quantities 
calculated  by  our  MFIE  computer  code  for  the  current  density  induced  by 
the  same  source  on  a flatly  capped  cylinder  having  the  same  length  and 
diameter  as  the  tube  of  the  experiment  are  presented  in  figures  6b,  7b, 
and  8b.  In  these  figures  we  employ  the  symbol  K rather  than  J_  to  denote 
the  current  density  in  order  to  conform  to  the  notation  of  the  measured 
data.  The  scales  of  our  calculations  were  adjusted  to  the  scales  of  the 
experimental  data  by  using  three  numbers,  a multiplicative  factor  for 
each  set  of  the  magnitude  comparisons  and  an  additive  factor  for  the  set 
of  phase  comparisons.  These  three  numbers  were  determined  by  forcing 
one  point  of  the  experimental  data  to  match  one  point  of  our  calculated 
data  on  only  one  curve  of  each  of  the  three  sets  of  curves.  The  reason 
we  include  figure  5 in  this  paper  is  to  show  that  there  is  only  a minimal 
measured  effect  of  capping  the  tube,  thus  justifying  our  comparing  our 
capped  tube  calculations  to  the  uncapped  measurements.  It  should  be  noted 
that  our  comparisons  with  the  data  were  for  h * 36  cm  while  data  in  figure  5 
corresponds  to  h ■ 84  cm.  It  is  possible  that  capping  the  shorter  tube 
could  have  a greater  effect  on  the  measured  surface  distributions.  This 
could  account  for  some  of  the  differences  between  the  experimental  data 
and  the  calculations;  however,  as  can  be  seen  the  difference  is  already 
quite  small. 

Now  it  is  necessary  for  us  to  discuss  the  frequency  at  which  we  made 
the  comparison  between  our  calculations  and  the  experimental  data.  The 
normalized  value  kh  • 1.5  it  determines  the  frequency.  First  we  note, 
without  scaling  h to  missile  or  aircraft  size  dimension,  that  the  comparison 
was  made  well  beyond  the  primary  resonance  of  the  cylinder.  Next  we 
mention  that  if  h is  taken  to  be  in  the  10  to  20  meter  range,  the  frequency 
scales  to  the  20  to  10  MHz  range.  In  this  regard,  we  claim  that  our  MFIE 
approach  can  only  perform  better  as  the  frequency  is  decreased  according 
to  the  explanation  contained  in  the  previous  section. 
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APPENDIX  A 


SYMMETRY  RELATIONSHIPS 


In  this  appendix  we  derive  certain  properties  of  the  surface  current 
density  induced  on  a perfectly  conducting  body  possessing  one  or  more  sym- 
metry planes.  The  symmetry  properties  are  true  when  the  body  is  illuminated 
by  a plane  electromagnetic  wave  of  a particular  direction  of  propagation  and 
polarization  as  we  will  explain  shortly.  These  properties  are  useful  in  that 
they  reduce  computation  time  and  also  give  insight  as  to  the  distribution  of 
the  induced  current  density  on  the  perfectly  conducting  body. ( See  reference  7 
for  a one  symmetry  plane  analysis  used  in  a different  context .) 

We  divide  this  appendix  into  two  parts.  In  the  first  part  we  consider 
a body  with  three  planes  of  symmetry  with  an  incident  wave  vector  k_  perpen- 
dicular to  a plane  of  symmetry  and  the  electric  field  parallel  to  an  axis 
of  symmetry.  We  derive  the  symmetry  relationships  for  a circular  cylinder, 
because  a circular  cylinder  is  relevant  to  this  report,  but  analogous 
properties  can  similarly  be  derived  for  any  body  with  three  planes  of  symmetry. 
These  relationships  involve  the  surface  current  densities  at  points  symmetric 
to  planes  of  symmetry  other  than  the  plane  perpendicular  to  the  k vector. 

In  the  second  part  we  consider  a body  with  only  one  plane  of  symmetry  with 
an  incident  wave  vector  perpendicular  to  the  plane  of  symmetry  and  the 
electric  field  parallel  to  a suitably  defined  axis.  We  show  that  the  sur- 
face current  densities  at  points  symmetric  to  the  plane  of  symmetry  are 
simply  related  to  each  other  as  the  frequency  oo  ■*  0 (magnetostatic  limit). 

Before  we  tackle  each  part  in  detail  we  present  certain  important 
results  whose  derivation  can  be  found  in. reference  3.  Starting  with  the 
magnetic  field  integral  equation  for  the  surface  current  density 

1 JJrJ  “ Imc^  + / dSo  <A-1) 

Js 

7.  Baum,  C.  E.,  Interaction  of  Electromagnetic  Fields  with  an  Object  Which 
has  an  Electromagnetic  Symmetry  Plane,  Interaction  Note  63,  Air  Force 
Weapons  Laboratory,  March  1971. 


37 


and  assuming  that  the  xy-plane  is  a plane  of  symmetry  for  the  body,  one  ob- 
tains the  following  set  of  integral  equations 


1 ,+  / + .. 
2 1 (E  ) 


-inc ('+)  + f ^+(^+;4)  ' J+(£o}  dS° 


S 

+ 


| J~(r+)  - J~nc(r+)  + J K_(r+;r^)  • j’<r*)  dSQ 


S 

+ 


(A-2) 


(A-3) 


where 


J+(r+)  + | [j(r+)  + Rz  • J(EZ  • r+)J 
J (r+)  = | £ J (r+)  - R 


R • J(R  • r+) 
z - -z  - J 


(A-4) 


(A- 5) 


r+,  r+  are  radius  vectors  at  points  with  z S 0,  R is  a reflection  operator 
- ~o  =z 

defined  by 


„ _ „ A A 

R = I - 2 e e 
-z  z z 


(A-6) 


I is  the  identity  operator,  S+  is  half  the  surface  defined  by  z > 0,  J?nc 
are  defined  through  the  incident  current  density  J ^ by  expressions  analt 
to  equations  (A-4),  (A-5)  and 


±.+  +»  . . + -K  _ , + + x 

K (r  ; r ) - A(r  ; r ) + B(r  ; r ) 
= - o - — o - — o 


A(r+;  r+)  “ K(r+;  r+) 
— “ o = o 


B(r+;  r+)  = K(r+;  R • r*)  • R 

— O — — 20  ^ 


(A-7 ) 
(A-8) 
(A-9) 
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Notice  that  r+  » x,y,z(iO),  jl  • r+  =*  x,y,-z. 


The  above  relationships  show  that  instead  of  solving  equation  (A-l) 
for  points  over  the  entire  cylindrical  surface,  we  can  solve  equations 
(A-2)  and  (A-3)  at  points  over  half  the  surface  and  calculate  the  current 
density  at  any  point  by  inverting  equations  (A-6)  and  (A-5),  i.e. 


J(r+)  = J+  + J~ 


J (R  . r+)  = R . (J+  - J”) 


With  the  above  as  background  information  we  now  consider  part  one  in  detail. 


1.  THREE  PLANES  OF  SYMMETRY 


The  geometry  is  depicted  in  figure  A1  where 

„ A -ikz 
H.  » - H e e 
-inc  o y 


(A- 10) 


A A , 

We  define  the  unit  surface  vectors  t,  s such  that 
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( A- 1 1 ) 


where  n is  the  unit  outward  normal.  Thus 
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(A- 13) 
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Jt(x,y,z)  » Jt(x,-y,z)  =*  Jc(-x,y,z)  - Jt(-x,-y,z)  (A-15) 

Jg(x,y,z)  » -Jg(x,-y,z)  « -Jg(-x,y,z)  = Jg(-x,-y,z)  (A-16) 

everywhere  on  the  cylinder  (walls  and  caps). 

These  relationships  give  us  a lot  of  information  concerning  the  distribu- 
tion of  the  surface  current  density.  For  example,  along  the  intersection  of 
the  y * 0 plane  with  the  cylinder,  J (x,o,z)  * 0.  Similarly,  along  the  inter- 
section with  the  x » 0 plane  the  azimuthal  component  Jg(o,y,z)  is  zero.  Notice, 
however,  that  equations  (A-15)  and  (A-16)  involve  points  symmetric  to  either 
the  x * 0 or  the  y * 0 planes  not  the  z = 0 plane.  Part  two  deals  with  symme- 
tries across  the  z = 0 plane  in  the  limit  u;  “ 0. 

First  we  consider  points  symmetric  to  the  y » 0 plane.  The  following 
preliminary  calculations  are  necessary. 

a.  Walls 


Let  us  calculate  the  incident  current  density  J c and  the  auxiliary  currents 


where  R^  • r » x,-y,z. 


With  the  aid  of  figure  A2  we  see  that 
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If  we  recall  definitions  (A-4)  and  (A-5)  we  see  that 


-inc(x,^,Z)  " Jinc(x*y’z) 


WX’y’Z>  - 0 


b.  x > 0 Cap 


With  the  aid  of  figure  10  we  see  that 
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.1  


and 


Therefore, 


R • J.  (R  . r)  - J.  (r) 
=y  -inc  =y  - -inc  - 


Jinc(x,y,z)  = -inc 


Jinc(x,y,z)  = 0 


( A- 2 1 ) 
( A-22 ) 


We  can  show  that  equations  (A-21)  and  (A-22)  are  also  true  for  the  x < 0 cap. 


Thus,  we  have  shown  that  at  every  point  where  equation  (A-3)  is  defined 
the  source  term  Jinc  is  zero.  (Notice  that  we  apply  equation  (A-3)  for  the 
y » 0 plane  of  symmetry  rather  than  the  z - 0 plane  for  which  it  was  originally 
presented . ) 


Equation  (A-2)  is  an  integral  equation  that  has  a unique  solution  and 
consequently 

J (x, y , z)  - 0 


( A- 2 3 ) 


everywhere  on  s+. 


From  equation  (A-5)  we  see  that 


J(r)  - R • J(R  • r) 

=y  ~ =y 

By  considering  the  s and  t components  of  equation  (A-24)  we  will  show  that 


(A-24) 


Jt(x,y,z)  - Jt(x,-y,z) 


J (x,y,z)  - -J  (x,-y,z) 

3 S 


( A- 2 5 ) 
(A-26) 
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everywhere  on  the  cylinder  (walls  and  caps). 
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c.  Walls 

(1)  t-Components 

If  we  take  the  inner  product  of  equation  (A-24)  with  t(r)  we  obtain 

J (r)  * t (r)  • R . J(x,-y,z)  (A-27) 

t - =y 

Recalling  that  t « e^  on  the  walls  we  see  that 

c(r)  • R - t(r)  - t (R  • r) 

-y  =y 

and  from  equation  (A-27) 

Jt(x.y.z)  * Jt(x,-y,z)  ( A- 2 5 ) 

(2)  s-Components 

From  equation  (A-24)  we  obtain 

J (r)  ■ s(r)  • R • J(x,-y,z)  (A-29) 

s - - =y 

With  the  aid  of  figure  A3  we  see  that 

S(r)  • |y  - [_sy (r)  + «z(r)]  , R, 

- -Jy( + 5z(r) 
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and  equation  (A-29)  gives 


J (x,y,z)  « -J  (x,-y,z) 
s s 


c.  Caps 


(1)  x > 0,  t-Components 


Writing  equation  (A-27)  for  this  case  and  with  the  aid  of  figure  A3 


we  obtain 
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and  equation  (A-27)  gives 


Jc(x,y,z)  - Jt(x,-y,z) 


(2)  x > 0,  s-Components 


Again  we  use  equation  (A-29)  and  figure  A3 
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and 


Jg(x,y,z)  » -Jg(x,-y,z)  (A-26) 

We  can  show  Chat  equations  (A-25)  and  (A-26)  are  also  true  for  the  bottom 
cap  x < 0. 

The  proof  for 

•Jc(x,y,z)  a Jt(-x,y,z)  (A-30) 

Jg(x,y,z)  - -Jg(-x,y,z)  ( A- 31) 

that  relate  the  current  density  components  at  points  symmetric  to  the  x * 0 

plane  is  similar  to  the  proof  we  gave  for  the  y = 0 plane.  In  the  x * 0 case 

we  can  show  that  j"*"  “0  everywhere  instead  of  J = 0 as  in  the  case  of 

inc  “inc 

y » 0 but  otherwise  the  proof  follows  the  same  lines.  Equations  (A-25), 

(A-26),  (A-30),  and  (A-31)  are  the  same  as  equations  (A-25)  and  (A-16). 


5.  INTEGRAL  EQUATIONS  WITH  FULL  REDUCTION 

We  will  now  employ  equations(A-15)  and  (A-16)  along  with  the  reduction 
scheme  given  by  equations  (A-2)  through  (A-9)  to  derive  a pair  of  integral 
equations  for  two  fictitious  current  densities,  defined  only  over  the  part 
of  the  surface  of  the  cylinder  that  corresponds  to  the  first  octant,  i.e. 
x 2 0,  y - 0,  z 2 0 and  show  that  from  a knowledge  of  the  current  density 
over  the  surface  of  the  cylinder  that  corresponds  to  only  two  octants  (x  ? 0, 
y ? 0,  z 2 0 and  x 2 0,  y 2 0,  z - 0)  we  can  calculate  the  current  density 
at  any  point  on  the  remaining  six  octants.  As  we  mentioned  earlier,  the  sur- 
face Integrals  in  equations  (A-2)  and  (A-3)  are  evaluated  over  S+  which  is 
defined  for  z 2 0.  By  applying  the  same  reduction  scheme  we  can  substitute 

++  f ■■ 

equation  (A-2)  by  two  integral  equations  for  J , J and  equation  (A-3)  by 
two  integral  equations  for  J +,  J by  considering  the  symmetries  about  the 


J 
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x • 0 plane  and  evaluating  the  resulting  surface  integrals  overs  defined 
by  x 2 0,  z i 0 . We  can  go  one  step  further  and  consider  symmetries  about 

'OO 

the  y - 0 plane.  Then  the  original  equation  (A-l)  for  J is  transformed  into 

f f f | | - |-  i — | .... 

eight  integral  equations  for  J_  , J_  , J_  , J_  , J ^ ^ , J_  , J_  , J_ 
with  kernels  K * * * » K^'  , KH  K*”  , K ' ' , K K +,  K respectively.  The 
surface  integrals  will  be  evaluated  over  s , , , defined  by  x 2 o,  y ? 0,  z > 0. 
By  repeated  application  of  defining  equations  like  equations  (A-4)  and  (A-5) 
we  can  derive  the  following  relationship 
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+ yR  • J(R  • r)  + aSR  • R • J(R  • R • r) 

»Z  — -.Z  — = x =sy  — = x =y  — 


+ ayRv  • _5Z  • J (JX  ’ —z  • L?  + SyRy  • J^2  • • J}2  • 

+ aByR  • R • R • J(R  • R ■ R • r)  I (A-32) 

=*x  =y  3z  — =x  *y  =z  — J 

where  on  the  left  hand  side  a, S.Y  are  + or  - and  on  the  right  hand  side  they 
are  equal  to  +1  or  -1.  Also  notice  that  in  deriving  equation  (a-32)  the  x - 0 
plane  symmetry  was  considered  first  followed  by  the  y * 0 plane  symmetry  and 
finally  the  z ■ 0 plane  symmetry  but  the  order  is  not  important. 


We  want  to  show  that  of  the  eight  currents  defined  by  equation  (A-32) 
six  are  identically  equal  to  zero.  We  will  do  so  by  considering  the  t and  s 
components  of  equation  (A-32).  By  following  a procedure  similar  to  the  one 
that  led  to  symmetry  relationships  (A-15)  and  (A-16)  one  can  show  that 


t<£)  • _§x  - * r) 

t(r)  • _Ry  - t(Ry  • _r ) 

^(r)  • R,  - t(R,  . r) 

— sz  — z — 


J 


( A-  3 3 ) 
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and  consequently 


t(r)  • R • R - -t  (R  • r)  • R ■ -'tCR  • R • r) 
— =x  =y  =x  — =y  =y  =x  — 


-t(R  • R • r) 
=x  =y  - 


since  R , R , R commute,  and 
=x  =y  =z 


Similarly 


and  consequently 


t (r)  • R • R - -^(R  • R • r) 

— =X  —Z  =X  =2  — 


t(r)  • R • R *t(R  • R • r) 
_ =y  =z  =y  =z  _ 


t (r)  • R • R • R » -t(R  • R • R ■ 
— —x  — y — z —x  — y — z 


A,  x 

s(r) 


A/  s 

s(r) 


s(r) 


R = s (R  • r) 
=x  =x  — 


R - -s(R  • r) 

-y  -y 


R * -s'(R  • r) 

= 2 =Z  _ 


(A-34) 


(A-35) 


(A-36) 


s(r) 

. R 

. R 

“ -s ( R 

• R 

• £> 

=x 

=y 

=x 

=y 

> 

s(r) 

• R 
=*x 

. R 
=>z 

- -'s(Rx 

• R 
= 2 

• I* 

► 

(A-37) 

s(r) 

. R 

. R 

A . _ 

- s(R  • 

R . 

r) 

~y 

25  Z 

=y 

> 

s(r) 

. R 

. R 

• R “ - 

■s(R  • 

R • R • r) 

* 

—x 

=y 

27 

=x 

=y  =*  - 
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Before  we  consider  the  t and  s components  of  equation  (A-32)  we  also  recall 
equations (A- 15)  and  (A-16)  which  we  can  rewrite  as 


J„(r)  - J(R  . r)  • J (R  • r)  » J (R  • R . r) 
c — c sex  — — c =y  — — c =x  =y  — 

JfR  . r)  - J (R  . R . r)  - J(R  . R . r)  - J (R  . R . R . r) 

C Z — t =aX  s2  — t assy  aZ  — t =X  =*y  3*2  — 

HA-38) 

J (r)  - - J (R  . r)  - -J  (R  . r)  - J (R  . R . r) 

s_  s=x—  s — y — — s =x  — y — 

j 

J (R  . r)  - -J  (R  .R  . r)  - -J  (R  . R • r)  * J (R  . R . R . r) 
s _z  — s =x  — z — s —z  — — s — x ^ —2  — 

In  view  of  equations  (A-33)  through  (A-38),  equation  (A-32)  gives 

■ j |q  + 8 -a  - aS)  Jt(r)  + (y  + 8y  - ay  - a8y)  J£(RZ  • r)j  (A-39) 
J*  “ 5 [d  + 8 - a - a6)  Jg(r)  - (y  + Sy  - ay  - a8y)  Jg(R2*  r)J  (A-40) 


From  equations  (A-39)  and  (A-AO)  we  can  see  that  only  J and  J are 
non-zero  and  that 


i[jt(r)  ♦ Jt(S2  • r>] 

I [-Vi’  - Viz  ' X>J 
I - Viz  ' x>] 
\ [vi> + Viz  ' I>] 


(A-41) 


Thus  we  can  evaluate  J (r)  at  any  r on  the  entire  cylindrical  surface  if  we 
-+*  ^ — 

know  J T Symmetry  relationships  (A-38)  show  that  this  is  possible  if  we 
know  J(r)  and  J(R2  • r)  where  r ■ (x  2 0,  y 2 0,  z 2 0) . Thus  from  a knowl- 
edge of  the  current  density  over  the  surface  of  the  cylinder  corresponding 
to  only  two  octants  (x  2 0,  y 8 0,  z 2 0 and  x > 0,  y 2 0,  z 1 0)  we  can  cal- 
culate the  currenc  density  at  any  point  on  the  remaining  six  octants.  From 
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equation  (A-41)  we  obtain 


d.  (r) 

- J-** 

+ d“‘ 

t — 

t 

t 

d (R 

• _£)  “ 

-++ 

J 

t —z 

t 

d (r) 

-++• 

- d 

+ d~" 

S — 

s 

s 

d (R 
s — z 

• r)  = 

-(J's 

(A-42) 


Equations  (A-42)  give  J(r)  and  J(R  . r)  (r»  (x  > 0,  y > 0,  z > 0))  in 

— + ==  -f  | 

terras  of  J . In  order  to  write  down  the  equations  satisfied  by  J 
and  J[  we  must  calculate  the  corresponding  source  terms  and  kernels. 

One  can  show  that  the  incident  current  density  also  satisfies  relationships 
(A-I5)  and  (A-16)  or  their  equivalent  ones  given  by  equations  (A-38).  Thus 


Jl£,t  = 1 [Jinc,t(^  + Jinc,t(4  *_r)] 


* i H s (r)  sin  kz 
o y - 


d.  “ 7 I d.  (r)  + J.  (r  r)l 
inc,s  2 L inc,s  — mc,s  — J 


H t (r)  cos  kz 
o y - 


Jinc,t  " 2 [Jinc,t-)  “ Jinc,t(|z  ‘ L5] 


-H  s (r ) cos  kz 
o y — 


J,"*-  “ 7 Td,  (r  ) - d . (R  • r)l 

inc,s  2 |_  inc,s—  inc,s=z 


* -i  H t (r  ) sin  kz 

o y*- 
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By  repeated  application  of  equations  (A-7)  through  (A-9)  one  obtains 


«I5Io>  ’ • I*>  • 5*  + • Io)  • £y 


± K(r  ;R  . r ) • R - K(r;R  . R . r ) • R . R 
= _ — z _0  =z  =.  _ =x  =y  _o  =x  =y 


? K ( r ; R . R • r ) . R . R ± K(r;R  • R • r )•  R • R 
— — «x  — z — o =*x  =z  =*  — =y  =z  — o =y  =z 


K(r ; R • R • R • r ) • R • R • R 
=»  _=x  =y  =z  _o  =x  =y  =z 


(A-43) 


-+± 


In  order  to  exhibit  the  integral  equation  for  £ in  component  form  we  must 
calculate  the  inner  products 


Recalling  that 


where 


we  have 


s(r) 

. K_++ 

A,  \ 

' «(Oi 

3(£) 

- 

— o 

A,  . 

— H- 

A , . 

A , v 

t(r) 

II?"! 

' 3(Io)* 

t(r) 

K(r;  Tq)  - n(r)  x [.VGCrjr^)  x i J. 


VGCr-.r^)  - Q(R)R 
ikR 

Q(R)  - e ^ (-1  + ikR) 

4rrR 


R - r - r , R - R 
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N(p, 4>o)  - x p)  • Rl  - Qo < ®2  * p)  ‘ R2  " ^3^3  x P>  * ^3 


and 


V$4  * + V*5  * P)  * *5  * V*6  * p)  ' *6  (A-47) 


- Q7<£7  x £)  • R7  + Qg  ($g  x $)  . Rg 
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(A-48) 


A A 

h ’ Vlo5 

$2  “ ix  * $1 

t,  - X • ♦, 

*4  * j,  • $i 
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5 -ly 
*8  -|, 
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R 
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=z 
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(A-49) 


The  explicit  forms  (A-46)  and  (A-47)  define  our  notation  for  M(p,$o)  and 

Mp»$0)  which  admittedly  is  not  impeccably  clear.  To  obtain  M(s,-sq),  for 

example,  we  substitute  p by  s(r)  and  $ (i  ■ 1,  2,  8)  by  -s^  where  s is 

the  unit  vector  s evaluated  at  r >*  (x  5 0,  y 2 0 , z 2 0) ! The  integrals 

—oo'oo 

are  evaluated  over  the  part  of  the  cylindrical  surface  defined  by  x - 0. 
yo*  0,  zQS  0. 

Thus  the  existence  of  three  symmetry  planes  and  broadside  incidence  (see 
figure  Al)  allows  one  to  reduce  the  problem  of  calculating  the  current  density 
on  the  cylindrical  surface  to  equations  (A-44)  and  (A-45)  which  involve  inte- 
grating over  one  eighth  of  the  cylindrical  surface  whereas  the  original  equa- 

— +’l 

tion  (A-l)  is  defined  over  the  entire  surface.  Once  we  know  J and  J 

we  can  employ  equation  (A-42)  to  evaluate  J(r)  and  J(Rz  • r)  where  r ■ (x  > 0, 
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y > 0,  z 2 0).  The  current  density  at  any  point  is  then  obtained  by  using 
symmetry  relationships  (A-38).  If  the  k vector  is  still  perpendicular  to 
the  xy-plane  but  the  electric  field  has  x and  y-components  we  can  consider 
the  two  components  separately  and  employ  the  relevant  symmetries  for  each 
case.  The  final  current  density  is  then  given  by  superposition  of  the 
current  densities  corresponding  to  the  two  electric  field  components.  If 
the  incident  wave  has  arbitrary  direction  and  polarization  symmetry  rela- 
tionships (A-38)  are  no  longer  valid.  One  can  still  reduce  equation  (A-I) 
to  integral  equations  over  the  x z 0,  y 2 0,  z > 0 part  of  the  surface  and 
solve  for  the  eight  current  densities  J"4.  From  equations  (A-32)  one  can 
solve  for  the  real  current  density  J evaluated  at  points  over  all  eight  oc- 
tants in  terms  of  the  J4±±  current  densities.  For  a specific  plane  wave 
some  of  the  symmetry  relationships  given  by  equation  (A-38)  may  hold  and  this 
would  reduce  the  number  of  non-zero  J±-±  current  densities. 

2.  ONE  PLANE  OF  SYMMETRY.  MAGNETOSTATIC  LIMIT 

Consider  a perfectly  conducting  body,  possessing  a plane  of  symmetry 
xy,  illuminated  by  a plane  wave  with  the  wave  vector  k perpendicular  to 
the  plane  of  symmetry  and  the  H vector  parallel  to  the  y-axis  as  depicted  in 
figure  A4.  The  x-axis  is  chosen  conveniently,  for  example  for  aircraft 
it  is  the  axis  of  the  fuselage. 

We  will  show  that  in  the  limit  uj  » 0,  where  to  is  the  radian  frequency, 
the  surface  current  density  at  points  symmetric  to  the  xy  plane,  are  simply 
related  to  each  other.  The  specific  relationship  will  be  presented  shortly. 

The  surface  components  at  the  surface  current  density  J will  be  denoted  Jg 
and  J^,  that  is  at  each  point  on  the  surface  we  have 

J - J s + J t (A- 50) 

- s c 
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A 

1 

I 


I 


A A 

where  s,  t are  orthonormal  surface  vectors.  In  order  to  construct  an  ortho- 
normal triad  n,s,t  where  n is  the  outward  unit  normal  to  the  surface,  we 
slice  the  body  with  a plane  perpendicular  to  the  x-axis  and  define  s as  the 
unit  vector  tangent  to  the  intersection  curve  with  the  body.  At  any  point 
P on  the  surface  the  unit  normal  n is  perpendicular  to  s because  s is  a sur- 
face vector.  The  ^ vector  is  then  defined  as  the  cross  product  of  n and  s, 
that  is 


A 

c 


A „ A 

n x s 


( A- 5 1 ) 


The  unit  vector  t is  a surface  vector  at  P because  it  lies  in  a plane  perpen- 
dicular to  n,  that  is  it  lies  in  the  tangent  plane  at  P.  (Figure  A4  shows  the 
components  of  the  triad  n.S,^  at  several  points  symmetric  to  the  xy-plane). 

Now  that  we  have  defined  s and  c we  can  present  the  relationships  between 
the  current  density  components.  These  are: 


J (r+)  = J (r  ) 

S — S — 

Jt(r+)  - -Jc(r  ) 


) 


J 


( A- 5 2 ) 


where  r and  r are  points  symmetric  with  respect  to  the  xy-plane,  that  is 


(x,y,z) 


r » (x,y,-z) 


z > 0 > 


J 


(A-53) 


j The  proof  of  relationships  (A-52)  involves  three  steps: 

(a)  Use  of  certain  symmetry  arguments  through  which  the  usual  magnetic  field 
Integral  equation  for  J(r)  is  substituted  by  two  integral  equations  for 
J*  and  J on  the  positive  (z  > 0)  half  of  the  surface  S where 

* 
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(A- 54) 


I [i<r+>  1 *, 


and  R is  a reflection  operator  about  the  symmetry  plane 


R » I - 2 e e 
z = z z 


( A-  5 5 ) 


(I  is  the  identify  operators.) 


(b)  Proof  that  at  u)  * 0 the  source  term  J^nc»  for  the  J+  integral  equa- 
tion, is  zero  everywhere  on  the  (z  2 0)  surface  of  the  body  and  con- 
sequently ,I+  * 0. 

(c)  Inner  multiplication  of  equation( A-54) f or  J+  with  s and  t and  use 
of  geometrical  properties  of  these  unit  vectors  to  finally  show  the 
validity  of  equation  (A-52). 


We  now  present  the  above  three  steps  in  detail.  The  magnetic  field 
integral  equation  for  J is 


j J(r) 


Ji(r)  + J K(r, 


Eo}  ' J-(t- o)dSo 


(A-l) 


where 


J.  (r)  - n(r)  x H (r) 
-inc  — — -inc  — 


( A-5b ) 


and  S is  the  surface  of  the  perfectly  conducting  body. 

Equation  (A-l)  as  we  mentioned  earlier  can  be  transformed  into  the  following 
pair  of  equations: 


1 \ 
2 J_  (O 


Jtnc(-+)  +/  K~(r-i  0 * J- (-ro}  dSc 


(A-2 ) 
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1 i'<£+>  ‘ Jtac<£+>  + / iV!  V ■ i 


dS„ 


(A-3) 


± + + 

where  K_  (r  , rQ)  are  defined  by  equations  (A-7)  through  (A-9). 

Next  we  calculate  J*  (r+)  from 
inc  — 


J+.  (r+>  - \ J,  (r+)  + R • J.  ( R . r+) 
inc  — 2 [—inc  — ~z  —inc  — z — J 


( A-57 ) 


Jinc(f+)  “ "(I+)  X”inc(I+> 


- a(r+)  x 3 1 H 

- y J e 


-ikz 


(A- 38) 


and 


where 


R . J . (r  ) 
_z  inc  _ 


Rz  .-|-[n(r')  x ej  Ho  e1^} 

[\(r~)  * *y  * Vr"’  ' %]  VT 
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n “ n + n + n 
-x  — y — z 


From  figure  A4  we  see  that 


n (r")  x e - n (r+)  x e 
—X  — y —X  — \ 


n (r  ) x £ - - n (r+)  x e 

-z  - y -z  - > 


. 
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t 

1 


and  from  Che  definition  of  Rz  given  by  equation  (A-55) 


. +.  A 

- n (r  ) x e 
_x  _ y 


R • f n (r+)  x e 1 » - n ( 

_z  [_  -x  _ y J -x 

R . f n (r+)  x £ - n ( r+ 

=z  [ -z  - y J -z  - 


) * e 


Using  these  relationships,  equation  (A-59)  can  be  written  as 


R • 'J.  (r  ) 
—z  — inc  — 


[n(r+)  x j Hoe 


and  from  equation  (A-57) 


’L<r+)  ■ - 1 [V>  - s ] <•' 

= i Hq  ^n(r+)  x ey  j sin 


-ikz  ikz. 
e - e ) 


As  k -*■  0,  equation  (A-60)  gives 


lim  J+  (r+)  = 0 everywhere  on  S 
_inc  — +• 


and  from  equation  (A-2) 


lim  “ 0 everywhere  on  S+ 


From  equations  (A-54)  and  (A-62)  we  then  see  that 


J(r  ) - -Rz  • J(r  ) 


which  is  the  vector  relationship  between  the  current  densities  at  r and  r 
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In  order  to  show  the  first  of  equation  (A-52)  we  form  the  inner  product  of 
equation  (A-63)  with  s(r+): 

s(r+)  * J(f+)  “ Js(i+)  “ -s(f+)  ' Rz  • J(r")  (A-64) 

The  right-hand  side  of  equation  (A-64)  can  be  simplified  by  expanding  s(r*).  R 

=z 

s(r+)  . _Rz  *|^sx(r+)  + sy(r+)  + sz(r+)j  . J*z 

- 0 + s (r+)  - s (r+)  ( A—  65) 

-y  -z 

From  figure  A4  we  see  that 

s (r+)  - -s  (r  ) 

-y  - -y  - 

s (r+)  - s (r~) 

_ Z __  Z _ 

and  equation  (A-65)  gives 

s(r+)  • Rz  - -s(r  ) (A-66) 

Thus  the  right  hand  side  of  equation  (A-64)  gives,  by  virtue  of  equation 
(A-66), 

Jg(r +)  - Jg(r")  Q.E.D. 

The  second  of  equations  (A-52)  can  be  shown  similarly  by  inner  multiplica- 
tion of  equation  (A-63)  with  *t(r+)  and  noticing  that  (see  figure  A4) 
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APPENDIX  B 

A PERFECTLY  CONDUCTING  ELLIPSOID  IN  A MAGNETOSTATIC  FIELD 


i 


I 


r 


In  this  appendix  we  calculate  the  surface  current  density  induced  on  a 
perfectly  conducting  ellipsoid  immersed  in  a magnetostatic-like  field,  i.e. 
we  assume  that  the  frequency  of  the  incident  electromagnetic  plane  wave  is 
high  enough  to  cause  negligible  penetration  but  low  enough  to  decouple  the 
magnetostatic  and  electrostatic  interactions  (see  reference  2 page  5 for  a dis- 
cussion). The  current  density  induced  due  to  the  electrostatic  interaction 
will  not  be  considered.  Thus  we  will  solve  the  magnetostatic  problem  depict- 
ed in  figure  B1  where  the  incident  magnetic  field  is 


H. 

— inc 


e 

y 


(B-l) 


and  the  normal  component  of  the  total  magnetic  field  on  the  surface  of 
the  ellipsoid  (b<  a<  c)  vanishes. 


First  we  briefly  explain  the  meaning  of  the  ellipsoidal  coordinates 
£,n.G»  defined  by  the  following  relationships. 


a2+S 


2 

b“+£ 


c2H 


> £ > -b 


( B—  2 ) 


x2  2 z2 

— 2 +2  + ~2  “ 1 

a +n  b +r,  c +n 


2 2 
-b  > t|  > ‘•a 


(3-3) 


y* — + -S—  + -4—  - 1 -a2  > c > -c2  (B-4) 

a +£  b +£  c +5 

The  first  family  of  surfaces  represents  confocal  ellipsoids  defined  by 
5 ■ const.  (5*0  corresponds  to  our  perfectly  conducting  ellipsoid.) 
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2 

The  second  (n  ” const.)  corresponds  to  hyperboloids  of  one  sheet  (b  + n < 0) 

2 2 

and  the  third  (£  » const.)  to  hyperboloids  of  two  sheets  (b  + £ < 0,  a + 

C < 0)  (fig.  B2)  all  confocal  to  the  ellipsoids.  Solving  equations  (B-2) , 
(B-3),  and  (B-4)  simultaneously  for  x,  y,  z we  obtain 


x 


y 


z 


+ f(£  + a2)  ( n + a2)  ( C+  a2)  1 

L (b2  - a2)  (c2  - a2)  J 

. I~(£  + b2)  (n  + b2)  (c  + b2)  1 
L (c2  - b2)  (a2  - b2)  J 


+ f(g  + c2)  (n  + c2)  (c  + c2)  1 
L (c2  - a2)  (c2  - b2)  J 


(B— 5) 


(B-6) 


( B— 7 ) 


In  order  to  understand  the  geometrical  significance  of  the  ellpisoidal 

coordinates  we  trace  how  the  above  conicoids  come  into  being.  From  equa- 

2 

tion  (B-2)  we  see  that  for  £ > -b  all  three  forms  are  positive  and  the 

resulting  surfaces  are  confocal  ellipsoids  ranging  from  a sphere  at  infinity 

2 2 2 2 2 2 
for  £ -*  » (x  + y + z -*■  £ ) to  an  elliptical  disk  with  semi-axes  a - b , 

2 2 2 2 
c - b lying  in  the  xz-plane  for  £ = -b  + 6 (6  0).  As  £ (which  we  call 

2 2 2 2 2 2 
H for  the  range  -a  , -b  ) passes  from  -b  +5  to  -b  - 6 the  sign  of  the 

second  term  in  equation  (B-3)  becomes  negative  and  the  resulting  surfaces 

2 2 

are  hyperboloids  of  one  sheet.  For  n » -b  - 5 (6  •+■  0)  the  hyperboloid 

degenerates  into  the  region  in  the  xz-plane  that  lies  outside  the  elliptical 
2 2 

disk.  For  n « -a  +5  (5  -►  0)  the  hyperboloid  is  flattened  into  the 

2 2 2 2 2 

region  in  the  yz-plane  "inside"  the  hyperbola  -y/ (a  -b)+z/(c  -a)-l. 

2 2 2 2 2 2 
As  n (which  we  call  £ for  the  range  -c  , -a  ) passes  from  -a  +5  to  -a  - 6 

the  first  two  terms  in  equation  (B-4)  become  negative  and  the  resulting 

2 2 

hyperboloids  now  have  two  sheets.  For  n ■ -a  - 6 (6  -*■  0)  the  correspond- 

ing hyperboloid  is  the  region  in  the  yz-plane  outside  the  hyperbola 
2 2 2 2 2 2 

-y  /(a  - b ) + z /(c  - a ) ■ 1,  i.e.,  it  has  two  separate  sheets. 

2 2 

Finally,  as  ( ■ -c  +6  (6  -*■  0)  the  two  sheets  are  flattened  into 

the  entire  xy-plane,  i.e.,  the  two  sheets  coalesce.  The  above 
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discussion  shows  thac  the  correspondence  between  the  ellipsoidal  coordinates 
and  spherical  coordinates  (see  figure  B3)  is 


! 

! 

J 


C ^ r 

n *-*  <j>  (b-8) 

p 

c «-*•  0 


In  this  appendix  we  are  primarily  interested  in  determining  the  current 

density  at  the  z ■ 0 intersection,  i.e.  at  points  along  the  arc  of  the  ellipse 

2222  2222 
x /a  + y /b  » 1.  If  we  set  z * 0 in  (B-3)  and  combine  it  with  x /a  + y /b  * 

l we  obtain 


2 2 2 2 2 
(x  + y )n  ■ n(a  + b ) + n 

2 2 

The  solution  q =*  0 does  not  lie  in  the  range  -a  ,-b  and  consequently 

2 2 2 2 

n - x + y -(a  +b)  (B-9) 


If  we  define  0 such  that  (see  figure  3) 
x = a cos  $ 
y ■ b sin 

we  can  derive  the  following  useful  relationship 

2 2 2 2 
n ■ -(a  sin  4>  + b cos  4>) 


0-10) 


(B-l  I) 


I 


which  we  will  use  later.  (We  can  also  obtain  equation 
2 

5 ■ -c  to  either  equation  (B-5)  or  equation  (B-6).) 
formulation  and  solution  of  the  problem.  The  incident 
be  derived  from  a scalar  potential 


( B— 1 1 ) by  setting 
Now  we  turn  to  the 
magnetic  field  can 
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<J)  « -H  y 

inc  o J 


The  induced  or  scattered  magnetic  field  can  also  be  derived  from  a scalar 
potential  satisfying  Laplace's  equation 


H =■  -V<t> 
—sc  sc 


Thus  the  total  magnetic  field  is  given  by 


H = -V($.  +4  ) * -V$ 

— inc  sc 


and  is  such  that 


H * -V4>  • n 
n 


at  C * 0 


(n  is  the  unit  vector  normal  to  the  surface  of  the  ellipsoid.) 


If  we  recall  equation  (B-6),  equation  ( B— 1 2 ) can  be  rewritten  as 

r.-  . l2w  . . . 2.  T h ^ 


H * b2)(n  * b2)  (r-  » bV 

o 2 2 2 2 

L (c  - b ; (a  - b ) J 


- Af(C  )f(n  )f(C  ) 


where 
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f(p)  ■ P + b‘ 


p - C.n 


( B— 17) 


A - ? 


E.2-*2,  <c2-b2,p 


In  view  of  equation  (B-16)  and  the  boundary  condition  (B-15)  the  scattered 
potential  <t  should  have  the  form 


d>  - Bg(Of(n)£(0 

sc 


where  the  functional  form  of  g(£)  will  be  determined  by  requiring  that  $ 


satisfy  Laplace's  equation. 


From  equations  (B-15),  (B-16)  and  (B-I8)  we  find 


B * -A 


df/dg 

dg/dC 


$ - Af(n)f(0<f(S)-g(0 


i"  dEVd^l  \ 
- o / 


The  two  surface  components  of  the  surface  current  density  J are 


• 5 * 0 


j . a !_  w 

n “c  3c 


, £ - 0 


where 
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II  ■»' 


1 [(n-C)  (n-ol  4 

2 R(n) 


£ - 0 


(B-23) 


where 


1 [(C-S)  (c-n)]  15 

2 R(?) 


C - 0 


(B-24) 


and 


[< 


R(p)  “ | ( p + a2)  (p  + b2)  (p  + c2) 


(B-25) 


Thus 


- f’(n)f(c)  [w(f,g)/  [dg/ac ] J 


- f(n)  f'(c)  [w(f,g)/  [dg/dc]] 


5 - 0 


( B— 26 ) 


f'(p)  - j — ^ . , p - Cfn 

(p  + b2)4 


(B-27) 


and  W(f,g)  is  the  Wronskian  of  f and  g: 


U(f.g)  - f - 8 *§ 


(B-28) 


Next  we  determine  g(£).  Laplace's  equation  in  ellipsoidal  coordinates  has 
the  form  (ref.  8,  page  59) 


(H  - O R 


£ aT 


(n») 


+ it  [ Rn  — 


3<t 


(Rn%) 


♦(e-n)tc  w [\ Tr)-° 


(B-28) 
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Substituting  $ in  equation  (B-28)  by  its  form  given  by  equation  C B— 18) 
sc 

we  obtain  the  following  equation 


(B-29 ) 


Noticing  that  f(£)  also  satifies  equation  (B-29)  we  invoke  a well-known 
result  for  a second-order  linear  equation  that  allows  one  to  obtain  a 
solution  if  an  independent  solution  of  the  same  equation  is  known,  i.e. 


g(0  - f( 


5)  f 2 d' 

J f (V  Rc 

(q  + bVf  :--j  -ii.nl,! 

J,  (£  + b ) \(C  + a 


(B-30) 


~2  2 
) (£  +b  ) (£  + c ) 


The  scattered  field  is  due  to  localized  currents,  i.e.,  $(£)  must  vanish  at 
£ ■ “.  This  is  secured  by  making  the  upper  integration  limit  in  equation 
(B-30)  infinite.  Now  we  are  in  a position  to  evaluate  fdg/d£j£  » 0 and 
the  Wronskian  at  £ - 0.  From  equation  (B-30) 

i&l  » _L_  ? d£  ^ _ T JL 

d£’k  * 0 2b  J0  (£+b2)  V (£+  a2)  (£  + b")  (£  + c")  ab2c 


1 


2ab  c 


2 (ao  ~ 2) 


( B—  3 1 ) 


8.  Stratton,  J.  A.,  Electromagnetic  Theory.  McGraw-Hill  Book  Company,  Inc. 
New  York  and  London,  1941. 
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where 


a “ abc 
o 


00 

f d- 5— 

Jn  (S  + b2)  V(S  + a2)  (5  + bZ)  a + 


7) 


(E 


The  Uronskian  can  be  evaluaced  by  recalling  that  g and  f satisfy  equation 
(B-29),  which  is  of  the  Sturra-Liouville  type,  i.e. 


W(f,g)  « C/R 


(B 


where  C is  a constant  to  be  determined.  If  we  evaluate  f and  g for  £•»<*> 
we  have 


f(0  - (C  + b2)1*  H. 


g(0  - f (C) 


oo 

I 


AL 


£(0  RC 


r'd  f — d_^_ 

’ J r5'2 

-a 


5/2  3 C 


and 


-4  , 2 : 1 


2 1 1 


-1/2 


w(f,g)  =.  ,r  -fj  > - f | 


-c  -3/2  - c r,/2 


Thus 


C - -1 


and 


(B 
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-32) 


i-33) 


-34) 


We  can  now  rewrite  equations  (B-26)  as 


1 f ( n+a 2 ) (n+c2)  (c+b2) 


[(c2-b2)  (a2-b2)]i 


n ' ao  - 2 [(c2 


( C.+a  ) ( :+c  ) (v+ b ) 


[(c2-b2)  (a2-b2)]^  L 5 


where  the  + sign  corresponds  to  y<  0 and  the  - sign  to  y>  0. 


Let’s  now  evaluate  the  current  component  at  points  z = 0,  x = a cos$  , 

2 2 2 2 2 
y « b sin  $.  In  this  case  £ * -c  and  from  ( B— 11)  n = - (a  sin  $ + b cos  0 ) 

Thus, 


2H  b | , | 

, . . o I COS  $ 1 

J - ( Z * 0)  - - __  O 'J  O 2 0 U 

o (a^sin  $ + b^cos  $ ) 2 


Jn  (z  = 0)  - 0 


Equation  (B-36)  shows  that  the  surface  current  density  at  the  z = 0 inter- 
section is  perpendicular  to  the  z-0  plane.  To  translate  into  J2  we 
recall  equation  (B-ll)  and  that  J,  * - ( 1 /h  ) 3$/3n  , * — ( 1 /ha> )3$/3<J>p 

( ■ -J ^ , at  i • 0).  (As  one  can  see  from  figure  3 the  angle  $ is  not  the 

usual  polar  angle  4 such  that  x • p($  ) cos  $ and  y = p ( 4>  ) sin  $ . How- 

^ P P P P P 

ever,  tan  $ «*  -r  tan  .5  and  3/3$,  3/3$  have  the  same  sign.)  Thus  equation 

op  p 

(B-36)  gives 


0 i $ s ff/2 


Jc  “ “Je  " Jz  > 0 


ff/2  s $ S it 


j,  - j6  - -J  > ° 


it  < <+>  < 3ir/ 2 
3 it/2  * <J>  < 2tt 


Jt  ' -J9  ' Jz  < 0 


JC  ' J0  " -Jz  < 0 


Thus  Che  sign  of  J _ is  Che  same  as  Che  sign  of  cost] 


Parameter  aQ  given  by  equacion  (B-32)  can  be  expressed  in  Cerms  of  an 
ellipcic  inCegral  of  Che  second  kind.  This  can  be  accomplished  by  making 
Che  subscicudon 


cos<{> 


x + b 


x + c 


which  transforms  equacion  (B-32)  into 

2abc 

3o  * 2 .2.3/2 

(c  - b ) 


Li 


tan“9  d0 


2 2 k 

■k  sin  9K 


(B— 3 7 ) 


where 


cosij)  ■ b/c 


2 c2  - a2 
c - b 


Using  reference  9 No.  782.03  we  obtain 


- ^£251.  J 1 - k^in2*  - m&J]  ( B— 38 ) 

c sinJ6  L 1 - k ^ 1 - k J 


where 


E(<J>,k)  - 


L" 


2 2 

k sin  9 K d9 


(B-39) 


9.  Dwight,  H.B.,  Tables  of  Integrals  and  Other  Mathematical  Data,  The  Mac- 
millan Company,  New  York,  Fourth  Edition,  1964. 
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is  the  elliptic  integral  of  the  second  kind.  Certain  limiting  cases  are  of 
interest . 


1 


1.  INFINITE  CYLINDER  OF  ELLIPTICAL  CROSS  SECTION 


I 


oo 


Thus  for  an  infinite  elliptical  cylinder 


(a  sin  0 +b  cos  0) 


1 


J<j,  “ o 


(B-42) 


wi  th 


— inc  ■ H°  ^ 


and  0 defined  in  figure  3. 


Notice  that  Jz (♦  *0)  =*  H0(l  + a/b) . When  a =*  b,  Jz  = 2H0  cos0  a well-known 
result.  It  is  interesting  to  note  that 

Jz(<f>«0,  aj‘b)/Jz(0-O.a=b)  - (l+a/b)/2  > 1 


In  order  to  see  how  much  the  infinite  cylinder  solution  differs  from  the 
ellipsoid  solution  (which  can  be  made  to  look  like  a finite  cylinder  of 
elliptic  cross  section  for  c much  larger  than  a),  we  cast  equation  (B-38) 
into  the  following  form 


a 

o 


2a 

a+b 


2ab 


i * .2.  2, is,  2 

( 1— b /c  ) (a 


— r r(*’ 

-bz)  L 


k)  - 


(l-b2/c2)i5  ] 


(B-43) 


2.  PROLATE  SPHEROID 

If  we  recall  definition  (B-32)  for  aQ  and  set  a**b  the  resulting  integral 
is  identical  to  the  one  obtained  by  Sancer  et  al  (ref.  1)  for  a prolate  spheroid 
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immersed  in  a magnetic  field.  To  verify  that  expression  (B-38)  is  also 

2 

correct  we  must  consider  the  limit  carefully  because  if  we  sec  a»b,  i.e.  k“=l, 
we  obtain  the  indeterminate  expression  0/0.  Thus  we  set 

a2=b2(l+62)  , a=b(l+!562+0((S4)) 

and  noting  that 


b . ,.  b2  ^ 

— = cos*  , (1  - — ) = sin* 

c 


1 - k 


2 , 
tan  * 


we  obtain 


(1  - k2sin20)!4  = (cos 2 0 + — ~~  sin20  )‘s 

tan*-* 


cos0  (1  + — sin20  + 0(54)) 
2tan“* 


a „ (l  + ^2+  0(54)) 

0 sin3,* 


f 


c-— j-  cos<>  (1  + Sj  ,52  + 0(64)) 

6“ 


tan  * 

62 


sin*  + 


2tan"' 


, , . 1 , 1 + sin* 

— (-  sin*  + — In  


l - sin* 


0(5 


]} 
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APPENDIX  C 


NUMERICAL  SOLUTION  FOR  THE  MAGNETIC 
FIELD  INTEGRAL  EQUATION  FOR  AIRCRAFT 

1.  MODEL 

Our  model  for  an  aircraft  is  depicted  in  figure  1 . The  fuselage  as  well 
as  the  rest  of  the  aircraft  components  are  modeled  as  elliptical  cylinders  of 
major  axis  2a^  and  minor  axis  21k  where  i denotes  the  component  or  body  under 
consideration. 

2.  ZONING 

Our  aircraft  is  symmetric  about  the  xz-plane  and  as  we  explained  in  Sec- 
tion I we  can  utilize  this  symmetry  to  transform  the  integral  equation  for 
the  current  density  J into  two  integral  equations  for  J"  defined  over  half  of 
the  aircraft  (y2  0).  Thus  we  will  only  zone  the  aircraft  for  y>0. 

The  half  airplane  for  y>0  consists  of  four  sections  or  bodies.  Body  1 
is  the  fuselage,  body  2 is  the  wing  corresponding  to  y>0,  body  3 is  the  hori- 
zontal stabilizer  (y2  0)  and  the  vertical  stabilizer  is  body  4.  Each  body 
is  treated  independently  in  this  appendix.  The  intersections  are  treated  in 

Appendix  D. 

Because  we  treat  the  bodies  independently,  we  define  the  orthonormal 
triad  n,  s,  t accordingly,  i.e.  s corresponds  to  the  azimuthal  direction  de- 

A 

fined  by  the  angle  <j>  (fig.  3),  n is  the  normal  to  the  surface  and  t » 
n x s.  Analytically, 

3 r / 3 r 

s “ a*/  3$  (c~ 


where  r is  the  radius  vector. 

Thus  at  each  point  on  the  surface  of  a body  we  have  a pair  of  orthonorroal 
surface  vectors  s and  c which  define  two  orthogonal  directions  and  for  this 
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reason  we  introduce  a convenient  surface  coordinate  system  (<p,Z)  where  <f> 
corresponds  to  the  s-direction  and  Z to  the  t-direction  on  the  walls  but 
the  radial  direction  on  the  caps.  Each  body  will  be  divided  into  4>~s  trips 
and  Z-s trips.  The  $-s trips  are  bounded  by  i - constant  lines  and  the 
Z-strlps  by  <{>  ■ constant  lines.  For  each  zone  a reference  point  is  taken 
and  it  represents  the  central  point  of  a zone.  The  coordinates  of  this 
central  point  are  given  in  Section  4 of  this  appendix.  The  meaning  of 
<J>  and  i can  be  clarified  by  considering  a specific  body,  say  the  fuse- 
lage. Any  point  on  the  surface  of  the  fuselage  can  be  described  in 
terms  of  three  parameters:  <p,  x and  r where  x ■ x,  y = a^  r sin<J>, 

z * -b^  r cos$.  On  the  caps  (x  * 0,  x =»  r ranges  between  0 and 

1 and  on  the  walls  r * 1.  Now  i can  be  defined  as  follows.  On  the 
walls,  that  is  for  0 < x < l let  i * x.  For  x * 0,  0 £ r £ 1 let  Z - Z^ 

1 + r.  For  the  rest  of  the  bodies,  which  only  have  one  endcap  at  x * length 

of  body,  l - length  + r for  the  endcap. 

We  are  now  in  a position  to  exhibit  our  zone  numbering  scheme.  As  we 
mentioned  earlier  we  have  four  bodies  which  we  have  numbered  from  1 to  4. 

Zone  no.  1 is  assigned  to  body  no.  1 (the  fuselage)  and  corresponds  to 
the  t-strip  defined  by  <p  a 0,  tp  ■ > 0 and  i - 0,  i » Z'  > 0.  The 

subsequent  zones  are  numbered  in  the  direction  of  increasing  Z until 

we  reach  Z ■ 2.^  + 2.  Then  we  go  back  to  the  t-strip  defined  by  $ - <J>^, 

4>  “ 4>2  > and  2*0,  Z ■ Z'  >0  and  the  subsequent  zones  are  numbered 
in  the  direction  of  increasing  Z until  we  reach  Z ■ 2^  + 2 and  so  forth. 

When  we  have  covered  the  fuselage  we  continue  with  the  wing  (body  no.  2) 
following  the  same  procedure,  i.e.,  using  the  same  numbering  scheme  in 
the  local  <p,Z  space.  Figure  Cl  illustrates  the  numbering  scheme  we  just 
outlined.  (Figure  Cl  should  not  be  interpreted  as  providing  any  informa- 
tion with  respect  to  the  relative  sizes  of  bodies  1 and  2 or  the  size, 
number  and  uniformity  of  zones). 

).  '1ATRIX  EQUATIONS  FOR  CURRENT  DENSITY 

Throughout  the  following  discussion  ct  denotes  a zone  number,  i denotes 
. : 'lumber,  J denotes  a strip  index  defined  by  two  $ - constant  boundaries, 
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icrinq  scheme  for  zones  on  the  aircraft  component 
conqxjnents  are  shown:  the  fuselage  and  the  wing 


and  k denotes  a scrip  index  defined  by  two  1 * constant  boundaries.  For 
example,  zone  28  in  figure  16  corresponds  to  a = 28 , i * 1 , j = 2 , k = 11. 

The  relationship  between  a and  i,  j,  k will  be  exhibited  separately  for  each 
body  in  section  A of  this  appendix.  We  wish  to  evaluate  the  s and  t-com- 
ponents  of  the  current  density  J at  the  center  of  each  zone.  Thus  the  de- 
sired quantities  are:  Jg(a),  Jt(a).  The  system  of  equations  for  finding 
these  quantities  is: 


" Jis<a>  “ H [A(“.“0>  Js(cto)  + W] 

rv  ^ J 


- Att  J,  (a)  = T'  C(a,a  ) J (a  ) + D(a,a  ) J (a  )| 
It  L oso  otol 


(C-2) 


where  the  source  terras  are 


^l(a) 

x Hi(a)J  = -t(a) 

• iii(ct)  I 

£n(a) 

x H^ooJ  = s(a) 

• H.(ct)  J 

(C-3) 


and  is  the  incident  magnetic  field  given  by 

H » H exp  (ik  • r) 

— i — o o “ 


(C-A) 


and 


k * -k  sin0  cos* 
ox  o o vo  1 


k - -k  sin0  sin*  V 

oy  o o o • 


k ” -k  cos9 
oz  oo 


(C-5) 


are  the  wave  vector  cartesian  components  (fig.  2) 


8A 


To  find  the  cartesian  components  of  we  have  to  define  the  polarization 
direction  of  the  electric  field  E^:  If  we  consider  the  axis  perpendicular 

to  k and  lying  in  the  _k,  z«axis  plane,  say  x",  then  forms  an  angle  ir-$ 
with  x"  (fig.  2).  (E^,  and  x"  lie  in  a plane  perpendicular  to  k.  We 

can  now  calculate  the  cartesian  components  of  by  making  three  successive 
rotations  as  follows.  We  assume  that  the  k,  E^,  system  initially 
coincides  with  the  -z,  -x,  y system  and  we  bring  it  to  its  final  position 
by  first  rotating  about  z by  $ (rotation  matrix  A)  then  about  y'  by  9 
(rotation  matrix  B)  and  finally  about  z"  by  -<t>  (rotation  matrix  C)  (fig.  C2) . 
Thus 


or 


H * H (cos*  cos6  sind>  - sinl)  cos<J>  ) 
oxo  oop  op 

H * H (sintfc  cos6  sinA  + cosi)  cos,*  _) 
oy  o oop  To  p 

H “ -H  sin9  sin<{) 
oz  oop 


(C-6) 


The  matrix  terms  in  equation  (C-2)  are  given  by  the  following  expressions. 


A(ot,a  ) - - .56  +|  R(a,a  ) • (c  <s  )Q(a,a  )dS' 

* o aa  I o a a o 

° W 
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4 


B(a,aQ)  = / R(a,aQ)  -(c^  ta  )Q(a,ao)dS 


s ( a ) 

o 


C(a,a  ) = - I 

° l 


S(a  ) 
o 


R(a,a  ) • (s  . )Q(a,a  )dS 
o ct  01q  o 


(C-7) 


D(a,a  ) 
o 


■56  - / R(a,a  ) • 

aa  I — o 

° J 


)Q<“.%)dS 

o 


3(a) 

o 


Q(a,a  ) * -1  + ik  R(a,a  ) exp  ik  R(a,a  ) /R3(a,a  ) 

o[_  o °JL°  °J  ° 


R(a,a  ) **  r -r 

— o —a  — a. 


R(a,a  ) - Ir  -r 
o '-a  -a 


(C-8) 


As  we  mentioned  earlier  our  aircraft  is  symmetric  about  the  xz-plane, 
and  we  can  utilize  this  symmetry  to  transform  ( C—  2 ) into  a pair  of  equations 
for  two  fictitious  current  densities  defined  over  only  half  the  surface  of  the 
aircraft,  i.e.  y2  0.  These  equations  are 

-47tJ“  (a)  - T |A-(a,a  ) J“(a  ) + B±(a,a  ) (a  )]  'I 
is  *—*  L oso  otoj  j 


(a)  * y.  C (a, a ) J (ct  ) + D (a, a )J  (a  )|  ! 

it  4 L oso  o t oj  J 


( C— 9 ) 


where 


87 


J (a)  » J (a)  + J+(a) 
s s s 

J (-a)  - J“(ct)  - J+(a) 

s s s 

. _ (C-LO) 

Jj.  (a)  - Jt(a)  + Jt(a) 

Jt(-a)  = J+(ot)  - Jt(a) 


Equations  (C-10)  show  that  from  a knowledge  of  J'  over  zones  corresponding  to 
y - 0 one  can  calculate  the  real  current  density  *J  over  the  entire  surface, 
i.e.  over  zones  corresponding  to  a and  -ct. 


The  source  terms  are  given  by 


where 


J~g(a)  = -t(a)  • H"(oO  e ^ 
- s(a)  • H~(a)  e +v^ 


H*(a)  - -isin * (H  e +H  e ) + cos*  H e 
l ox  x oz  z y oy  y 

H . ( a)  - -isin  * H e + cos  * (H  e +H  e ) 
1 oy  y v ox  x oz  z 


( C- 1 1 ) 


( C— 12) 


9 ■ k z(a)  cose 
o o 


^ “ k x (a) sin 0 cos*  r 

o o Yo 

^ “ k^yCaJsinQ^siniJi^  > 

The  above  equations  are  derived  by  first  recalling  that 

if  (a)  » j * |y-  JjC-o)] 


( C— 13) 


(C-14) 


i 
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where 


-k  • r(±a) 


J.(±a)  * n(±a)  * H (±a)e 
— i — o 


If  we  define 


k • £(a)  ■ -(9+iJ^-t) 


and  use  the  following  relationships 


tx(ot)  ' Cx("a) 


cy(a)  ” -ty(_a) 


t (a)  = t (-a) 
z z 


sx C a.)  “ -s  (-a)  ' 


sy(a)  = sy(-a) 


s (a)  * -s, (-a) 
z z 


which  can  easily  be  demonstrated,  one  can  show  that 


- t » "lx'-’’] 

* [■..«*  1 

* [^ly(o>  1 V-^] 


and  equations  (C-ll),  (C- 12)  can  then  easily  be  verified. 


Finally,  we  give  the  expressions  for  the  matrix  terms  on  equation  (C-9) 


A-(a,ao)  - -0 . 56aao  + J [R+(a,aQ)  • (ta  x s^)  Q(R+) 

S(aQ) 

R“(a,ao)  • (to  x Ry  • sao)  Q(R-)]  dS 

B±(a,ao)  - J [R+(a,ao)  • (ta  x tM)  Q(R+)  ± r" (o,aQ) 

S(aQ) 

• (ta  X Ry  • tao)  Q(R")J  dS 

<T(a,ao)  m ~ J [R+(a,ao)  • (sa  x sao)  Q(R+)  ± R"(a,aQ) 
S(ao) 

• <VV  sao)  QdO]  ds 

D!(a,oc0)  - -0.56^  - j [R+(a,ao)  • <Sa  x Q(R+) 

S(«o) 

± R‘(a,ao)  • (3a  x Ry  • t^)  Q(R-)]  dS 

Q(R“)  * [-1+ik  R"(o,a  )]  exp[ik  R±(a,a  )]/[R±(a,a  )]3 
o o o o o 

R+(a,a  ) • R(a,a  ) - r(a)  - r(a  ) 
o o — — o 

R (a, a ) - r(a)  - R • r(a  ) 

— o — -y  — o 


(C-15) 


(016) 
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* 
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In  Section  5 of  this  appendix  we  present  a detailed  calculation  for  A“, 

+ + + 

B~ , C“  and  D“. 

4.  DEFINITION  OF  COORDINATES  OF  CENTERS  AND  BOUNDARIES  OF  ZONES 

In  this  section  we  give  the  coordinates  of  the  centers  of  zones  that 

are  not  adjacent  to  the  intersections.  Assuming  that  we  have  N.  trans- 

)ti 

verse  strips  (defined  by  Z m constant  boundaries)  on  the  walls  of  the 
ith  body,  Nri  transverse  strips  (defined  by  l - constant  boundaries)  in 
the  end  cap  of  the  ith  body,  Nsi  longitudinal  strips  (defined  by  0 ■ 
constant  boundaries)  on  the  ith  body  and  N zones  on  the  ith  body  we  can 
now  present  the  defining  relationships  for  the  coordinates  of  the  centers 
of  zones  that  are  not  adjacent  to  intersections  (Appendix  D) . 

a.  Fuselage  (Body  No.  1) 

One  can  easily  show  that 


a - a(i«l,j,k)  - (j-l)(N.  + 2N_  ) + 

*i  i 


(C-17) 


where  a,  i,  j and  k were  defined  in  Section  3 of  this  appendix.  The  index 

j varies  from  1 to  Ng^  and  k from  1 to  Nj^  + 2Nr^  depending  on  whether 

we  are  on  the  walls  (k  * l,...,No  ),  on  the  front  cap  (k  » Nn  + l,...,No 

1 *1  1 
+ Nr^)  or  the  back  end  cap  (k  - Nj^  + Nr^ Nj^  + 2Nr^). 

je[i,NSl],  kftl.N^] 


vf  M)  %-£(*-*) 


x(a)  - x^ 


s (a)  » 0 
x 


tx(a)  ‘ 1 


b.  cos0 

y(a)  - bj  altupj  ,yfa)  - 

a1  sin<|> 

z(a)  - -ajcos^  sz(a)  - 


t (a)  - 0 

y 


tz(a)  - 0 


(C-18a) 
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where 


N^(<J>j)  * (b2  cos2^  + a2  sin2$  )*^2 
je[l,NSl],  k - NAl  € [l,Nr  ] 

1 J S1 
j ) * (a^  sin2^  + b2  cos 2^) 1/2 


x(a)  ■ 0 


s (a)  =•  0 
x 


tx(a)  * 0 


b,  cos<J>, 

y(o.)  - h rfc  sln+j  .yta)  - 


a sin<|> 


z<«)  - -J  vk  cos^j  »z<a)  - -g-jp 


a sin<t>. 


V<c0  ■ "vv 


-b  cos4>. 

t!‘“) 


j6[l,N  ],  k - N - Nr  e [1,N  ] 

31  *1  1 rl 

\ ’ [r  (k  - \ - \ - 1)]1/2.  *rir(>-  i) 


W ■ (bl  C0SVj  + al  sin2^)1/2 


x(a)  - l. 


8x(a)  * 0 


b i cos4>, 

y(a)  - bt  rk  .1^  .y(«)  . 


a.  sin$> 


«<a)  - -at  rk  co.^  .,(<.)  - -^-y 


t (a)  - 0 
x 


-a  sin<J>. 
t (a)  “ — J 

V ' N^) 


b.  C03(p 

Cz(a>  " Nx  (<)>.) 


(C—  20) 
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b. 


Wing  (Body  No.  2) 


where 


I 


a - a(i-2,j,k)  - (j-l)  (N5  + N ) + k + w (C-21) 

*2  r2  1 

N^  is  the  total  number  of  zones  on  the  fuselage  (body  no.  1),  i.e., 

’ 2Nr  > Ns  ' 

1 S1 

je[l,N  ],  kf[l,N  1 
S2  *2 

(j  ~ i)  N2^j)  " (a2  + b2  cos^j^ 

■ sf  (k  - l)  + bl 
X2 


x(a) 

- xQ2  + a2  cos4»j 

sx(a) 

-a2  sin<p^ 

n2(4>j) 

t (a)  * 0 

X 

y(a) 

* yk 

sy(a) 

- 0 

ty  (°)  - 1 

z(a) 

■ -b2  sin<()j 

sz(°t) 

b2  cos<}>^ 

n2  (^  ) 

t (a)  - 0 
z 

(C— 22) 

jeti.H  1. 
S2 

k - 

N,  f 
2 

tl.N  ] 
r2 

♦j  “ 

360  / 

-1) 

H2(*j)  - 

f 2 .2.  , . 2 2.  .1/2 
(a2  sin  + b2  cos 

r,  ■ 

JL/k 

- N 

-I) 

1/2 

k 

[Nr  A 

2. 

*2 

2/ 

x(a) 

y(a) 

z (a) 

c. 

where  N 


x(a) 

y(a) 

2(a) 


X02  + a2rk  cos*j  sx(a) 


a2  sintf^ 


b„  cos0 , 

cx(a>  - - 


l2  + al 


sy(a)  - 0 


ty(a)  « 0 


"b2rk  sin*j 


b,  cos4> 


i 


S,(0l)  M (l  \ 

2 N2  (<J>j  ) 2 


Horizontal  Stabilizer  (Body  No.  3) 


t (a) 


a2  sin<j>. 

W 


(C-23) 


a - a(i*3, j ,k)  - (j-l)(N  + N ) + k + N + N,  (C-24) 

*3  r2  12 

N2  are  the  total  number  of  zones  on  bodies  1 and  2,  respectively. 


],  ke[l,N.  ] 
3 *3 


fM) 


N3 (<J>j ) = (b2  cos2<}>.  + a2  sin^j)1^2 


sf  (k-i)+*i 


x03  + a3  C0S*i 


s (a) 
x ' 


-a^  sinij)^ 
n3(4>j  ) 


tx(a)  - 0 


s (a) 

y 


ty(a)  - 1 


-b^  sin<}ij 


s2(a) 


b2  cos^ 

N3OM 


tz(a)  - 0 


(C-25) 
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J€[1,N  ].  k - N,  e [l.N  ) 

SL  j r3 


vf(-l) 


“ (b ^ cos2^j  + a3  sin2^^)1^2 


57  (k  - \ ■ i) 


1/2 


*(“)  * xQ3  + a3rk  cos0j  sx(a) 


sin<j> . 
N3(d)j) 


b_  cos<J> 

tx(a)  * " N3(*.)j 


y(a)  * + a1 


s (a)  * 0 

y 


ty(a)  = 0 


z(a)  * -b3rk  sln4>^ 


b,  cos<J) 

s (a)  — ■ .J 

z 


d.  Vertical  Stabilizer  (Body  No.  4) 


a.  sin$. 
t (a)  - -J 1 


N3(4>j) 


(C-26) 


a - ct(i-4,j,k)  - (j-l)(N  + N ) + k + N.  + N + N,  (C-27) 

^4  r4  1 L J 

where  N^,  N^,  N3  are  the  total  number  of  zones  on  bodies  1,  2,  and  3, 
respectively. 


je[l,Ng  ],  kffl.N^  ] 
4 4 


V 'iif  (3  'i)  ^4 ^ j 3 


,2  ,2.  , , 2 2.  .1/2 
(a^  sin  + b^  cos  4>j) 


4MH 


(C— 28) 
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x(a)  " XQ4  + a4  cos<<,j 


y(a)  - 8in<f>j 
z(a)  - zfc 


s (a) 
x 

sy  («) 

s2(a) 


*4  8in<<)1 

■4»,) 


b.  cos$, 

-i _1 

Vty 


t (a)  - 0 
x 

t (a)  - 0 

y 

tz(a)  - 1 (C-28) 


3€tl,N  ),  k - M 6 [1,N  ] 

3 4 *4  r4 


v**H) 


Vty 


(a^  sin2(J>j  + b2  cos2*))^)1^2 


■k  (k  ■ \ - 


1/2 


a,  sin$. 

x(0)  ’ x04  + a4rk  C°8*j  sx(0t) N^) 


tx(a) 


b.  cos0. 

_ J> 1 

Vty 


y(a)  - b4rfc  sin<() 


b.  cos$, 
8y(a)  “ Vty 


ty(a) 


a.  sin<f>. 

_ _i 1 


Vty 


z(a)  ■ + ax  sz(o)  - 0 t2(a)  - 0 (C-29) 

Equations  (C-17)  through  (C-29)  give  the  coordinates  at  the  centers 
of  zones  that  will  be  used  in  our  numerical  solution.  For  the  calculation 
of  matrix  elements  given  in  the  next  section  we  need  to  know  the  boundaries 
of  the  zones  over  which  we  integrate. 

For  zones  on  the  walls  we  can,  in  a straightforward  manner,  use  the 
corresponding  coordinates  for  the  centers  to  calculate  the  limits  of 
integration  for  the  matrix  elements.  Thus  for  the  fuselage  we  refer  to 
Xj^  and  <£j  in  equation  (C-18)fl  and  for  the  (i«l,j,k)  zone 
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4>J  " *2  ($2_  + 4>2) 


(018). 


and  similarly  for  the  other  components  of  our  model. 

For  zones  on  the  caps  the  0 ■ constant  boundaries  are  easily  obtained. 
Thus  for  the  front  end  cap  on  the  fuselage  we  refer  to  <J>^  in  equation  (C-19). 
and  we  can  readily  obtain  4^  - (180/Ng^)  (j-1) , <J>2  - (180/NSl)j.  For  the 
r ■ constant  boundaries,  r^  at  the  center  of  a zone  is  not  the  average 
of  the  r's  of  the  boundaries.  Instead: 


r~  (k  - \ - 0 ' 

*2  1 'J 

t (k  ‘ \)1  W2 


2 1.2,  2. 
rk  ’ 2 (rl  + r2> 


* _ 180  . . . _ 180  . 

h ~ (j  1}  *2  ~ j 


»k  - 2(*1  + *2>  (C-19>b 


and  similarly  on  the  back  end  cap  of  the  fuselage  or  the  caps  of  the 
other  components  of  the  aircraft. 

Before  we  go  on  with  the  calculation  of  matrix  elements  in  the  next 
section  of  this  appendix,  we  would  like  to  state  an  important  feature 
built  into  our  computer  code.  Equations  (017)  through  (029)  dictate 
the  rule  that  a given  zone  obeys  in  relationship  with  a neighboring  zone 
on  the  same  surface  (wall  or  cap).  To  allow  for  sufficient  nonuniformity 
for  experimenting  with  the  zoning  Che  code  provides  that  equations 
(017)  through  (029)  can  be  applied  sectlonally;  that  is  we  first  divide 
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the  walls  or  caps  in  sections  by  drawing  selected  <f>,  t,  r-boundaries 
and  then  within  the  sections  we  apply  the  rules  that  govern  equations 
(C-17)  through  (C-29).  For  example  we  can  divide  the  wall  of  the 
fuselage  in  three  sections  in  the  x-direction  and  two  sections  in  the 
^-direction  and  choose  any  number  of  zones  within  each  section  to  allow 
for  the  desired  nonuniformity  of  our  zoning  scheme. 

5.  CALCULATION  OF  MATRIX  ELEMENTS 

+ + £ + 

In  this  section  we  calculate  the  matrix  elements  A”,  C and  D“ 
in  equation  (C-15).  First  we  will  outline  the  method  of  derivation  and 
then  perform  a sufficient  number  of  calculations  that  will  allow  the 
reader  to  understand  how  equations  (C-62)  through  (C-92),  which  give  the 
matrix  elements,  were  obtained.  These  equations  are  still  valid  for  the 
self-terms,  i.e.,  when  a is  the  central  point  of  the  zone  over  which  we 
integrate  (see  Appendix  D for  zones  adjacent  to  instersectlons) . 

As  we  can  see  from  equation  (C-15)  the  calculation  of  matrix  elements 

+ £ + + + + 

A-,  B , C and  D"  involves  integrands  of  the  form  R • (p  x qQ)Q(R  ) and 

R * (f>  x R • 3 )Q(R  ) where  3 and  3 are  unit  surface  vectors.  The 
" o o 

integrated  variable  has  the  subscript  zero  and  runs  over  a particular  zone 
whereas  the  free  variable  corresponds  to  the  center  of  a zone  anywhere 
on  the  four  bodies.  All  four  bodies  are  elliptical  cylinders  and  as  we 
explained  earlier  the  surface  unit  vectors  s and  £ are  chosen  to  conform 
to  the  geometry  of  the  body,  i.e.,  3 is  defined  in  the  azimuthal  direction 

Ak 

defined  by  the  angle  $ (fig.  3)  and  t is  equal  to  n x s where  n is  the 
outward  unit  normal  to  the  surface.  To  facilitate  our  subsequent  calcula- 
tions, to  each  body  we  attach  a cartesian  coordinate  system  x^,  x^ 
such  that  x ^ ■ x^,  X2  ■ ar  cos$,  x^  ■ br  sin<J>  where  r - 1 on  the  walls 
of  the  body  and  0 £ r <_  1 on  the  endcaps.  The  correspondence  to  the 
global  coordinate  systems  xyz  is 


Fuselage 

Wing 

Horizontal 

Stabilizer 

Vertical 

Stabilizer 

X1  " x 

xi 

xi  - y 

*i  * 1 ] 

N 

1 

• 

N 

K 

x2  “ x * x02 

X2  " X “ x03 

x2  - x - xQ4  . 

\ (C-3Ca) 

x3  • y 

x3--z 

X3  ■ -z 

x3  - y J 

The  eadcaps  of  Che  fuselage  are  Chen  deCermined  by  ■ 0 and  x^  * of 
Che  wing  by  ■ a^  + &2»  of  the  horlzoncal  stabilizer  by  » a^  + £3 
and  of  Che  vercical  scabilizer  by  ^ - bx  + £^.  The  relacionship  beCween 
Che  angle  $ depicted  in  figure  1 and  $ just  defined  is 


Fuselage 

Wing 

Horizontal 

Stabilizer 

Vertical 

Stabilizer 

<j>  *■*  ♦ 

<p  p 

<p  <— ► <f> 

( P (j) 

(C-30b) 

As  we  mentioned  earlier  the  free  variable  in  Che  integrands  corresponds 

Co  Che  center  of  a zone  anywhere  on  Che  surface  of  a body  with  local  unit 
A ^ 

veccors  s^  and  tQ  . For  a number  of  subsequent  calculations  we  do  noc 

have  to  specify  the  form  for  8 or  t and  we  will  denote  them  by  the 

a a 

symbol  v.  In  order  to  calculate  the  matrix  elements  we  must  calculate 
Integrals  of  the  form 


I+  - / R+  • (*  x 3 ) Q(R+)  dS . 
J.  — o o 


r - / i"  • (*  x r • 


y ° 


q J Q(R~>  ds. 


(C-31) 


where  8 is  either  5 3 or  t 2 t . 

o ao  o ao  o 


We  start  with  I+. 


For  a zone  on  a body 


and 


R+(0  x 3)  - ^*1(v2s30  - v,s20) 

) 

i On  walls  and  caps 

(C-32a) 

‘ vl(Ax2830  - Ax3820) 

) 

R+(v  x t)  « Ax2v3  - Ax3v2 

| On  walls 

R+(0  x t>  - AXl(v2t30  - v3t2Q) 

) 

/ On  caps 

(C-32b) 

- v1(Ax2t3Q  - Ax3t2Q) 

1 

^1  ’ xla  * X1 

**2  ’ x2o  * x2  " x2a  ' ar  c08*’ 
Ax3  " x3a  ~ x3  “ x3a  " br  8in*» 


®20  “ “■  ain$/N($) 
s30  - b eoa<p/N($) 

N($)  - (a2  sin2<|>  + b2  cos2<f>)1/2 


/ 0 £ r 1 for  caps 
r - 1 for  walls 

On  walls  and  caps 


t2Q  » b cos$/N($) 
t3Q  • a sin<i/N($) 
t2Q  • -b  cos$/N($) 
t3Q  - -a  sin*/N($) 


| On  • x 


0 endcap  of  fuselage 
All  other  endcaps  (C-33) 


Notice  that  we  have  simplified  the  notation  and  - (x3,x2,x3)  instead 
of  (xio,x20,x30^  to  which  we  return  later. 

Next  we  need  the  differential  surface  element  dS  for  a zone  on  the 

o 

walls  or  the  caps  of  a body.  On  the  walls  dS  - dx.ds  where  ds  is  the 

O 1 

arc  length  in  the  ^-direction 
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i 


ds  - [(dx2)2  + (dx3)2]  / - (a2  sin2<|>  + b2  cos2<i>)1/2  d4> 


Thus 


dSQ  - N(iJ>)dx^d<J> 


on  the  walls 


(C-34) 


To  calculate  the  area  dSQ  on  the  caps  we  recall  that  x2  ■ ar  cos$. 


x^  ■ br  sin$  and  ■ r(a  cos<t>  e2  + b sin$  e^)  where  r^  is  the  radius 
vector.  The  area  dS  is  thus  given  by 


dS  * | ds  x dr  | ■ 

n 1 — o ' 


3r  3r 

—n 


•o  ~o 
3d>  X 3r 


d((>dr 


dSQ  ■ abrdrd$ 


on  the  caps  (C— 35) 


Let  us  now  evaluate  integral  I (s)  given  by  equation  (C-31)  for  a 
zone  on  the  walls  au4  the  caps  but  not  adjacent  to  an  intersection  (see 
Appendix  D) . 


a.  Zone  on  Walls 

We  rewrite  equation  (C-32a)  as 


R • (0  x s) • (xla  - x3)  MjCifr)  + ^(40 


(C-36) 


Thus 


I+  (s) 
w 


/ 


R+  • (v  x 3)  Q(R+)  dSQ  - I*x  + I+2  (C-37) 


where  w - wall  and 


Cl  - j M1(^)  N(<»  d4>  J 


^■2  X.  ■ X-  jv  n 

N($)  d4>  | (ik  R - l)e  *°R  dx.  (C-38) 

jo  l 


4xR 


‘11 


I2  - j M^*)  N<4»)  d$J 


12  ik  R “ 1 a p 
N($)  d4>  | ■ °-  y ■ e1RoR  dx1 


(C-39) 


4ttR 


‘11 
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If  we  notice  that 


Xla  ~ X1 


3R 

3x, 


and 


/ ik  R\ 


ik  R 
o 


ik  R ik  ik  R 
e ° + ~p—  e 0 - *_ 
R R2 


(ik  R - 1) 
o 


we  can  rewrite  equation  (039)  as 


**!<♦>  m) 

* i 


r ikoR2  ikoRl 
e e 


R-, 


R, 


d<p 


R(*l>  ' [Ka  ' xl)2  + (x2a  - x2)2  + (x3a  ‘ x3)2] 

Rj_  m R(x2  m 1 


1/2 


R2  " R^xl  " x12^ 


(C-40) 


(C-41) 


Integral  Jw2  does  not  lend  itself  to  such  a simple  treatment  as  1*^.  As  a 
first  step  we  observe  that 


dxl  _ X1  ~ xla 
R3  a2R 


where  R2  • a2  + (Xj^  - x1(j)2  and 


f elk°R  (X1  ’ xl>  ° fik  »(*,  * x,  )2  e 0 

(v  2 2 2 

x2  “ x2a)  " R ” a w*  can  rearrange  this  equation  to  obtain 


(C-42) 


ik  R 


ik  R 
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(C-43) 


The  Integral  In  the  angular  bracket  Is  not  elementary.  One  can  expand 
exp(ikQR)  in  a Taylor  series  and  evaluate  the  resulting  elementary  Integrals 
but  too  many  terms  will  be  needed  for  a zone  remote  from  the  zone  over 
which  we  integrate.  This  can  be  remedied  by  multiplying  and  dividing  the 
Integral  by  exp(ikQRo)  where 


R(x  - x^)  + R(x  - x12) 


(C-44) 


Thus 


f*2 

C2  * 4ttJ  M2m  N(*>  d* 


' ikQR2  ikQR1 

K(4>)  - — pT  (xla  ' x12)  rT  (X1o  " Xll) 
a z.  1 


iko  ikoR0rXl2  iko(R-R0) „ 

~T e J e dxi 


'll 


Now  we  can  expand  exp[iko(R-RQ) ] ! 


/Xl2  ik  (R-R  ) JA  fX12  [ik(R-Rn)]n 

6 °dvLj  dxl 

Xil  n“°  X11 


(C-45) 


(C-46) 
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To  determine  nm  we  observe  that  max|R  - Rq|  % “ xi2^2  and  reca^1 

that  we  require  on  the  order  of  ten  zones  per  wavelength.  This  translates 
into  k max|R  - R0|  being  of  the  order  of  unity  and  consequently  we  only 
need  few  terms  to  secure  sufficient  accuracy. 

If  we  use  the  binomial  expansion  for  (R  - RQ)n  we  obtain  integrals  of 
the  form  /Rn  dXj^  which  are  elementary.  They  can  be  evaluated  with  the 
aid  of  the  recursion  formula 


(C-47) 


b.  Zone  on  Caps 

We  rewrite  equation  (C-32a)  as 


R+-  (v  x s)  - AXl(v2s30  - v3s2Q)  - Vl[(x2a  - ar  cos*) (b  cos*) 

- (x3(J  - br  sin*)(-a  sin*)]/N(*)  - ^xl/v2s30  ” V3S20^ 

“ vl(x2aS30  - x3aS20)  + Vab/N(<),) 

R+  • (v  x s)  ■ K^(*)  + rK2(*> 


(C-48) 


and 


l+a>  - J R+  (V  x s)  Q(R+)  dS0  - l+cl  + I^2 


(C-49) 


c ■ n f2  %'»> 


r2  ikoR 

re-  - (ik  R - 1)  dr 
n 3 0 


(C-50) 


r* 2 r 2 

J - 


r2  2 ikoR 


(ikQR  - 1)  dr 


(C-51) 
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and 


R2  ■ (xla  - x^2  + (x2a  - ra  cos<t>)2  + (x3a  - rb  sin<}>)2 

- , \2  . 22 
f * (xla  - xl>  + x2a  + x3a 


g - ~2(^2a  3 cos^  + x3a  b sini^ 


2 2 2 2 
h ■ a cos  0 + b sin  $ 


In  order  to  simplify  equation  (C-50)  we  observe  that 


rdr 

n3 


-d  2gr-  + 4,£  = -dP(r) 

(4fh  - g ;R 


/ 


re 


ik  R 
o 


ik  R 


dr  - -P(R)e 


/' 


° + 1 HR)lk0elk°R  i(. 


Performing  the  algebra  in  the  integrand  on  the  right-hand 
(C-53)  we  obtain 


P 


ik  R 


o ik  R 

dr  - -P(r)e  0 + 


/(>?K 


ik  R 

e ° dr 


where 


A ■ 4fh  - g . 


■ f + gr  + hr2 
(C-52) 


; + 2hr)dr 

(C-53) 
side  of  equation 

(C-54) 

(C-55) 


105 


If  we  rearrange  equation  (054)  we  obtain 


•r2  r(ik  R - 1)  eik°R  ik 

~ dr  - P(r)  e 0 


ikoR  r2  2ikog  fr2  ikoR 

T J e dr 


h 

t+  . ab  f 

cl  47 T I *1 

J A 


(<P)  d« 


I 2gr  + 


..  ik  R 2 
o 

— e 


(056) 


The  last  integral  in  the  angular  bracket  can  be  evaluated  by  the  procedure 

outlined  in  connection  with  equations  (043)  through  (046).  I+.  can  also 

c2 

be  simplified  by  noting  that 


^ dr  - d [-rP(r)  + 


where  P(r)  is  given  by  equation  (053)  and  A by  equation  (055).  Thus 


r 2 ik  R 
r e 


dr  " _rp(r)+^A  e ° +/r(^  + ^) 


ik  R 

ik  e 0 dr 
o 


. ik  R , C o 

ir  + 8)  n o . . 1 I e . 

"hA  iko  e dr  + h R dr 


This  relationship  allows  us  to  write: 
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2r(gr  + 2f ) 2gR 

ik  R 
0 
e 

R h 

A 

r2 

rl 


ik  R 
o 


dr  + 


tko8 

hA 


ik  R 
o 


e 


(C-57) 


The  integrals  in  the  angular  brackets  can  be  evaluated  as  in  equation  (C-56). 

To  calculate  I+(t)  given  by  equation  (031)  we  follow  the  procedure 
employed  for  I+(s)  and  arrive  at  equations  (040) , (045) , (056)  and  (057) 
where 


M^O)  - 0 

M^)  - &x2v3  ~ Ax3v2 

K1(<’)  “ Ax1(v2C30  “ v3t20)  “ vl(x2aC30  " x3at20) 

(<t>)  * ±v^(a2  - b2)  sin<|>  cos$  (058) 


and  the  + sign  refers  to  the  x^  - x ■ 0 endcap  of  the  fuselage  and  the  - 
sign  to  all  other  endcaps. 

To  calculate  I~(t)  and  I (3)  given  by  equation  (C-31)  we  observe  that 
gy  refers  to  the  global  coordinate  system  and  its  relationship  to  the 
local  coordinate  systems  is 


Fuselage 


Wins 


Horizontal 

Stabilizer 


Vertical 

Stabilizer 


£y-a*: 


Sy-Sxl 


R 

-y 


R i 

— xl 


R - R , 
■7  — x3 
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Fuselage 


—y  qo  ■— xl  % * ^qol,qo2* 


R « r - 


4c2  ’ — o * (xlo  ‘ x10)  el  + (x2ct  ” x205  e2 


+ + x30)  ®3 


(C-59) 


Wing  and  Horizontal  Stabilizer 


*y  * % " Kl  * % " (‘<»ol^o2*qo3) 


* " • *xl  * lo  * (xla  + x10)  ai  + (x2a  ' x20)  h 


+ (x3a  " x30)  e3 


(C-60) 


Vertical  Stabilizer 


Sy  ‘ ' 2,3  ’ ' <,ol',,o2,"',o3) 


S •■£<.-  Sx2  - io  " <xla  ‘ >10>  81  + <*2o  ‘ x20>  *2 


+ (x3a  + x30>  S3 


(C-61) 


We  can  use  the  previous  relationships  to  calculate  I (t)  and  l"(s)  on  the 
vails  and  endcaps  of  the  various  bodies  by  following  the  procedure  employed 
for  the  calculation  of  I+(t)  and  I (s).  Instead  of  exhibiting  these  results 
we  present  the  final  expressions  for  the  matrix  terms  A-,  B“,  C~  and  D± 
defined  by  equation  (C-15) 
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■I— 


■ - 1 V 


[<■ 


liJ*t,8o*  + I2ij(t,so)J 


± [Ilij(t,So)  + I2ij(t,so)] 


B <a*°o>  " [Ilij(t,to)  + I2iJ  ^t,Co^j  * [Ilij(t,t0)  + I2ij(t,to)] 

C <o’o0>  * ‘ |[IliJ*8,So)  + I2ij(s,so)]  1 [Ilij(s’so)  + I2ij(s,so)j| 

“‘(“•V  ' - I ««o  - { + I21j(s>to)] 


± [Ilij(s,to)  + 


3| 


(C-62) 


4*  4* 

Before  we  present  the  defining  equations  for  1”^  and  1“^  an  explanation 

of  the  notation  is  in  order.  Index  i runs  from  1 to  4 and  refers  to  the 

four  bodies.  Index  j ref  to  either  the  wall  (w)  or  the  endcap  (c)  of 

a body.  In  the  fuselage,  j refers  to  both  endcaps  (and  also  the  wall). 

+ 

The  I~  in  equation  (C-62)  are  simplified  expressions  of  the  Integrals  in 
equation  (C-15)  which  are  evaluated  over  a zone  on  the  surface  of  a body 

characterized  by  (<t>lt<J>2),  (xioi»xio2^  on  the  walls  or  ($!>$<>)»  ^ri>r2^  on 
a cap  and  by  the  surface  unit  vectors  t and  sq  given  by  equation  (C-33). 
The  free  variable  a corresponds  to  the  center  of  a zone  characterized  by 
its  coordinates  xlo,  x2a>  x^  and  unit  vectors  t , sa(s]L,s2,s3) . 

The  transformations  from  the  local  coordinate  systems  x^,  x2,  x^  to  the 
global  coordinate  system  x,y,z  are  given  by  equation  (C-30). 
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A~  Matrix  Term.  Fuselage 


± b^  cosi p + tj  sinO 
M-(0)  - tx[±  a^  + x2a  bx  cosO-  #1  sino] 


iT($) 


lkoR2 


L “I 


(xla  - x102) 


.ikb*l 


“i 


± (xla  * X101} 


ik„  ikX  rxi 02  ik  »*-»*) 

e dx 


+ eikoRo  r 102 

*i  X101 


10 


R *x10*  " [(xia  “ X10)2  + (x2a"  *1  coa^Z  + \ •*■*#> 

*1  " R±(x10  * x101> 

*2  - »*(*!„  - «102) 

. 'o  ■ i<*i + <4> 


(C-63) 


(C-64) 


jji/2 


(C-65) 
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Ill 


hi 

sin2*  + b 

2 2. 
^ cos  * 

)l" 

f + 

g±r  + hr2 

(xla 

-L  + : 

2 

*3(1 

^x^  b^  sin* 

- 2x 
* 2a 

a^  cos* 

2 2 2 
cos  * + b^ 

sin2* 

4fh 

-<»V 

(i-2 

I12w(t,so) 


4i  / “!<♦> 


W'-V 


'♦l 

.♦ 

'* 


ik  R*  ik  r:  ** 
e 0 2 e ° 1 


■£  J 


d* 


r*2 

J *£<♦>  K±C<P)  d* 


M{(*)  - t2b2  cos*  + t3  »2  sin*,  (M+  - M“) 

M*(*)  - [«2b2  - (X2a  b2  cos*  + X3a  a2  sin*)]tlf  (>£  - M^) 


(C-70) 


(C-71) 


(C-72) 
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K“(*) 


" a2  _±  (xla  T x102)  * _±  (xi«  ± xim 


la  " "101 


ik  iknr  /*"io2  tk  (ir-io 

+ — e I e dx 

a9 


10 


101 


R-(*10>  ■ f * xm>  + <x-v,  - ao  cos*)2  + (x 


la  T "10 


2a 


3a 


*1  " R (x10  " x101) 


R2  " R (x10  “ x102) 


± i + + 

R0  - -|<RI  + R2> 


A~  Matrix  Term.  Wing  (i-2.1-c) 


+ a2b2  C 2 + + 

Ii2c(t*so)  • ~trj  Kl(v  HI<*> d* 

♦l 

+ a2b2  C 2 + + 

I22c(t,so)  * ~ 4ir  I K^>  H2W  d<t> 

\ 

*!<*)  • [(xlQt  * xc)(C2  b2  C08<^  + C3  a2  sin*) 

* ci(x2a  b2  co*<>  + ^ a2  *in*)]/N(*) 

KI(*)  - + t.  a,b,/N(*),  (K+  . k7) 


1 

sin*)2] 
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x ■ b,  + l 


l + (endcap  of  vertical  stabilizer) 


N($)  ■ (a2  sin2$  + b4  cos2$)*^2 

4-  4*  ? 

R ■ f + g“r  + hr* 


f " (xla  * xc)2  + x2o  + x3 


8'  - -2x2a  a4  cost  * 2x3a  b4  sin$ 
h - a2  cos2$  + b2  sin2$ 


±*2 


A*  - 4fh  - (g~) 


B Matrix  Term.  Fuselage  (i-1.  1-w) 


Iilw(t'to)  " 0 


I21w(t*to> 


l^2  * 
-r- 1 m; 
4ir  I 2 


m:<4>)  ka(*>  dt 


M2(^)  *[t3(x2a“  ax  cost)  - c2(x3a  + bl  sin<^>]N(d>) 
K~ (t)  are  given  by  equations  (C-65). 

4- 

B~  Matrix  Terns.  Fuselage  (i-1.  1-c) 


*lbl  f*2  t t 

• J v*> d ♦ 


I21c(t,to) 


Vlf*3 

4*  J 


k2(«)  h2(^>  d t 


0 


(C 


(C 


:-84) 


:-85) 


:-86) 


:-87) 
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K^($)  ■ (-1)"  [(xla  “ xcX±t2  a^  sinij)  - t3  b.^  cos4>) 
-tl  ("  x2a  al  sin<*»“  x3a  bx  cos«)]/N(W 


K2(0)  «±(-l)n  t^  (a^  - b^j  sin$  cosd>/N C4>) 


(C-88) 


3 ^ 

H^($)  and  H2($)  are  given  by  equation  (C-68)  through  (C-70) , n »0 
for  the  x£  ■ 0 endcap  and  n » 1 for  the  xc  ■ 2^  endcap. 


B~  Matrix  Term.  Wing  (i-2.  1-w) 


I22w(t,to) 


f*2 

- «*(♦)  d<J> 


(C-89) 


M^(4>)  -lt3(x2o-a2  cos<)>)  - t2(X3a  - 


2(x3a  “ b2  sin4>>]  NC0>,  (M2  - M2)  (C-90) 


K “($)  is  given  by  equation  (073). 


B~  Matrix  Term.  Wing  (i-2,  1-c) 


I12c(t,to) 


a2b2  C2  t t 

- J «£(♦>  H'W)  d* 


I22c(t,to) 


J 1$(*)  H2(4>)  d<|> 


(091) 
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1 


4 

K^(<1>)  - [(xla  + xcHt3  b2  cos<j>  - t2  a2  sin<(>) 

“ ^(x^  b2  cos<<i  ~ x2a  a2  sin0)]/N(<j») 

* -t1(a2  - b2)  costp  simp,  (K+  m K~)  (O 92) 


+ 4 

H"(ij>),  H~ (4>) , x£  are  given  by  equations  (076)  and  (077). 

4 

B~  Matrix  Term.  Horizontal  Stabilizer  (i-^ 

Integrals  and  I^Ct.t,,)  are 

given  by  formulas  identical  to  those  for  the  wing  provided  that  a 2 is 
changed  to  a3>  b2  to  b3  and  x3  - y,  x2  - x - xQ3>  x3  - -z. 

4 

B~  Matrix  Terms.  Vertical  Stabilizer  (j»4.j»w) 


I"  (t,t  ) - 0 
14w  o 


(C-93) 


M^(4>)  - [(x2o  - a4  cost}) ) t3  - (x3a  f b4  sin<f>)  t2]  N(<{>)  (094) 

4 

K“(4>)  is  given  by  equation  (080). 

B~  Matrix  Terms.  Vertical  Stabilizer  (i»4,1»c) 


K*(4>)  H*(<fr)  d<J> 


k2(<j>)  H2(<|>)  d<p 


(095) 
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K^(4>)  * l(xict  “ x£)  (+*2  a4  sln<*>  + c3  b4  cos^) 

+ tl(±x2a  a4  sin<,)  " x3a  b4  cos<^)  ] /N(*) 

* ±t^(b^  - a^)  cos<|>  s in4> / N ( 4> ) (C-96) 

+ ■+■ 

H~C4>)  > xc  are  given  by  equations  (C-83)  and  (C-84). 

+ 

C Matrix  Terms 

Integrals  *2ij  are  identical  t0  those  calculated  for 

the  A terms  provided  that  tffl(m  = 1,2,3)  are  substituted  by  s^Cm  =*  1,2,3). 
+ 

D Matrix  Terms 

+ i 

Integrals  1“  (s,t  ),  1“  (s,t  ) are  identical  to  those  calculated 
+ lij  o Zij  o 

for  the  B_  forms  provided  tm(m  • 1,2,3)  are  substituted  by  sm(m  = 1,2,3). 
6.  SELF-ZONE  INTERACTION 

When  the  observation  point  r_  is  the  central  point  of  a zone  over 
which  we  integrate  we  talk  about  self-zone  interaction  or  self-zone 
terms.  The  integration  variable  r^  ranges  over  the  entire  zone  and 
consequently  at  the  central  point.  What  about  the  resulting 

singularities  in  the  integrand?  If  we  examine  the  A± , B~ , C*  and  D± 
matrix  terms  we  see  that  the  only  terms  of  concern  are  those  of  the  form 


(The  A , B , C and  D matrix  element  cannot  involve  self-terms  because 

R * r J*  .£_•)  Tbe  integration  is  over  a zone  on  the  cap  K+(tJ>)  has  the 

2 2 ^ 
form  ab/N(<p)  or  t^(a  - b ) sin<j>  cos$/N(<£).  But  when  r^  is  on  the 

cap,  = 0 and  the  potentially  troublesome  integral  is  identically 

equal  to  zero.  In  conclusion,  all  equations  giving  the  A*,  B±,  C~ 

D terms  are  valid  for  both  zone  to  zone  or  self-zone  interactions. 
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7.  INTERPOLATION 


In  Chis  section  we  explain  the  interpolation  scheme  that  allows  us  to 
evaluate  the  current  density  at  points  other  than  the  centers  of  the  zones 
on  the  surface  of  the  aircraft. 


As  we  discussed  in  Section  3 of  this  appendix  one  uses  equations  (C-9) 

to  evaluate  the  fictitious  current  components  J~  at  the  centers  of  the 

zones.  Notice,  however  that  the  variable  a in  equations  (C-9)  is  not 

restricted  to  represent  the  center  of  a zone.  Such  restriction  is  true 

+ 

only  for  Thus  we  can  use  equations  (C-9)  to  evaluate  Jg  at  any 

point  (see  section  9 for  points  near  edges  and  junctions)  on  the  surface 

+ 

of  the  aircraft  in  terms  of  J . at  the  centers  of  the  zones.  The  real 

s,t 

current  density  components  J can  then  be  calculated  with  the  aid  of 

s , t 

equation  (C-10) . 

8.  CHARGE  DENSITY 


The  charge  density  a can  be  calculated  with  the  aid  of  the  continuity 
equation 


iuxJ  - 7 • J 

s — 


(C-99) 


where  Vg  • is  the  surface  divergence  of  J^.  We  approximate  directional 
derivatives  by  central  difference  quotients.  We  can  therefore  write  the 
divergence  of  J as 


V 

s 


* i = 


3Jg(r)  3Jt(r) 
3s  + 3t 


J« (♦+*♦. t)  - J (<M<M)  J^.t+At)  - J (<fr,t-At) 

“ S , L L 

(AC)s  + (AC) t 


where  (AC)g  and  (AC)t  are  the  chords  between  points  (<(>+A<(>,  t) , (<J>— A4> , t)  and 

points  ($,t+At),  ($,t-At)  respectively.  In  the  above  equation  the  values 

of  J and  J are  calculated  through  the  interpolation  scheme  as  we  explained 
3 C 
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in  the  previous  section.  If  J and  J could  be  calculated  with  unlimited 

t s 

accuracy  and  it  were  possible  to  calculate  the  desired  differences  with 
similar  accuracy,  then  as  At  and  A$  approach  zero  the  error  of  the  approxi- 
mation to  7g  • J would  approach  zero.  However,  since  the  interpolation 
formula  is  an  approximation  to  equations  (6)  based  on  a finite  number  of 

current  density  values  we  understand  that  AJ  /(AC)  and  AJ  /(AC)  depend  on  the 

s s t t 

accuracy  of  J (<J>+A4>,t),  J (<t>-A<f>,t) , J ($,t4At)  and  J„(<|>,  t-At) . If  we 

s S t C 

introduce  the  relative  errors 


J®  - J 
s s 


at  “ 


Jt-  Jt 


where  the  superscript  e stands  for  exact,  then: 


numerical  divergence  = "V  • J"  » ....  . 

s - (AC)  (AC)  (AC) 


AJ®  AJ®  A(agJ®)  A(atJ®) 


s 


(AC), 


In  the  limit  (AC)  -*■  0,  (AC)  + 0 we  have 
s c 


AJ®  AJ®  5J®  9J® 

+ ..  - - — — + — — ■ io) O® 


(AC)  (AC)  5s  5s 

S b 


f J 


and 


9J  9a 

iuo  - "V  • J"  - i(ja®  - a -r-5-  - J®  -r-5-  - a 
s — s 9s  s 9s  t 


9J®  9a„ 

— - J®  — 
9t  Jt  9t 


9J  9a  9a 

m>m  ' 1“’  + IT  - “,>  + IT  + IT 
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J 


For  regions  where  the  error  varies  slowly  with  position,  J®(3as/3s)  + 
Jt(3at/3t)  can  be  smaller  than  the  first  term  in  the  square  brackets. 
Equation  (C-100)  then  shows  that  In  the  low  frequency  limit  the  numerical 
calculation  of  the  charge  density  can  be  very  inaccurate.  When  (AC)g  and 
(AC)t  are  finite  the  error  ag  will  be  given  by  equation  (C-100)  plus 
correction  terms,  but  the  low-frequency  problem  will  persist  even  though 
it  may  be  lessened. 

The  above  treatment  is  by  tw  means  complete  but  it  provides  insight 
into  the  difficulties  that  beset  the  numerical  calculation  of  the  charge 
density  with  the  aid  of  the  continuity  equation  as  required  in  the  work 
statement  for  this  proposal. 

Due  to  time  constraints  we  did  not  explore  a different  method  that 
is  based  on  solving  an  integral  equation  for  a whose  source  term  involves 
the  normal  component  of  the  incident  electric  field  and  an  integral  over 
J^.  The  equation  is 

<*o> 


foM 


’r)  ■ iwe  n(_r) 


<i>  - f[< 


k0G(R)  n(r) 


Using  the  same  zoning  scheme  as  for  the  solution  for  we  can  calculate  a 
at  the  centers  of  the  zones  and  subsequently  at  any  point  through  our 
interpolation  scheme. 
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9.  EDGE  AND  JUNCTION  BEHAVIOR 

« 

In  this  section  we  will  discuss  the  behavior  of  the  surface  currents 
paying  particular  attention  to  the  effects  of  discontinuities  in  surface 
normals.  We  define  an  edge  (junction)  to  have  more  (less)  air  than 
metal  in  the  vicinity  of  the  discontinuity.  We  begin  our  discussion  by 
rewriting  equation  (1)  in  component  form 


1 JS<E>  ■ + / W “e 

S 

+ /kb(^)  Jt< 


.(r  ) dS, 


2 Jt(£>  * Jit(i>  + / KC(£*o>  Js(^o>  dSc 

•'s 

+ f* d(X5Io)  Jt< 


.(r  ) dS. 


where  the  kernels  ka>  Kg,  Rc,  KQ  are  defined  as 


KA  - (r  - r ) • (t  x s ) F(k0R)/R3 


kb  ■ (£  - rQ)  • (t  x tQ)  F(k0R)/R3 
Kc  • -<r  - rQ)  • (s  x sq)  F(k0R)/R 


K0  * ’ Iq)  * (s  x tQ)  F(k0R)/RJ 


ik  R 

F(k0R)  - (-1  + ikoR)  ^4tF~ 


R " li-iol 


(C-lOla) 


(C-lOlb) 


(C-102) 


(C-103) 
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We  proceed  by  examining  the  behavior  of  the  kernels  and  their 
integrals  I„(r)  ■ /K(r;r  ) dS  . We  notice  that  each  of  the  kernels  are 

K.  — O 0 

unbounded  as  R approaches  zero,  and,  since  |f|  <_  2 the  kernels  are  well 

behaved  everywhere  else.  In  addition  I„(r)  will  be  bounded  whenever  the 

1+e 

triple  product  involved  decreases  as  fast  as  R , e > 0,  as  r^  approaches 
£.  If  the  surface  has  continuous  curvature  in  the  vicinity  of  the  reference 
point,  jr,  the  three  vectors  smoothly  approach  being  coplanar  making  e j>  1. 
Thus  our  kernels  vary  no  faster  than  1/R  and  consequently  the  kernels  have 
an  lntegrable  singularity  (see  reference  10  for  a rigorous  treatment  of  the 
case  that  results  in  e ■ 1) . However,  if  _r  is  near  a discontinuity  in 
curvature  the  above  reasoning  may  fail  and  I (r)  may  be  unbounded.  We  will 
now  specifically  treat  the  behavior  of  the  kernels  in  the  vicinity  of  a 
discontinuity  in  the  direction  of  the  normal  (i.e.,  for  edges  and  junctions) 
and  then  discuss  the  theoretical  and  numerical  consequences. 

We  present  our  arguments  by  considering  two  plates  Intersecting  at 

right  angles.  More  complicated  geometries  as  they  appear  in  our  model 

would  not  change  our  results  but  they  would  make  our  arguments  harder  to 

follow.  If  we  expand  F(kQ|ir  - r^|)  in  a MacLaurln  series  we  see  that  only 

the  zeroth  order  term  will  contribute  to  the  singularity  of  our  integrands. 

We  therefore  confine  our  attention  to  k *0. 

o 

Our  model  is  depicted  in  figure  C3.  We  choose  s on  both  plates  to  be 
parallel  to  the  Intersection  of  the  plates  and  let  t « n x s.  If  we 
treat  our  intersection  as  a junction  our  surface  vectors  become 

81  • 32  • V S1  “ V a2  " ‘V  *1  " -*x  and  *2  " -gz  (c"104) 

We  place  the  reference  point  on  plate  1 and  limit  our  attention  to  a 
rectangular  patch  (denoted  by  Sj)  on  plate  2 which  Includes  the  point  on 


10.  Marin,  L.  and  R.  W.  Latham,  Analytical  Properties  of  the  Field 
Scattered  by  a Perfectly  Conducting.  Finite  Body.  Interaction 
Note  92,  Air  Force  Weapons  Laboratory,  January  1972. 
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Figure  C3: 

! 


Two  metallic  plates  intersecting 
at  right  angles  and  description 
of  n,  s,  t. 
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plate  2 closest  to  the  reference  point  since  only  that  patch  can  cause 
a singular  Integral.  The  Integrals  over  that  patch  are  then 


.s  „b  ~ 


K.  dS  = A 
A o 


-iff  <C1  * *2>  —12  , , 

* W J J ^3 dy: 


.3  3 


Kb  dSo  = B 


-1  C C (tl  x t2)  ' -12  . 

SiJ  J JS ^2^2 


s b 


Kr  dS  = C 
u o 


1 C C (S1  x S2>  * ^12  . . 

J i? ^ 


3 Jj 


Kn  dS  = D 
D o 


(s1  X t2) 


dz2  dy2 


(C-105) 


where  — 2*  and  we  ^ave  chosen  y^  ■ 0. 

s ^ equals  §2  making  C Identically  zero  and  t^  x t2  equals  - making 
Kg  antisymmetric  in  y2;  thus  both  B and  C are  well  behaved.  For  the  geom- 
etries of  our  model  it  can  be  shown  that  both  Kg  and  Kc  vary  no  faster  than 
1/R  as  jr  ■+•  even  in  the  presence  of  a discontinuity  and  consequently 
are  well  behaved. 

The  Integral  D is  also  well  behaved  but  for  a much  more  subtle  reason. 
Carrying  out  the  integrations,  we  have 


This  expression  is  bounded  for  all  values  of  the  parameters.  The  saving 
feature  of  this  integral  is  that  the  proportionality  factor  is  not  an 
integration  variable.  This  feature  will  also  be  present  for  more  compli- 
cated geometries. 

Unlike  the  situation  for  D,  the  proportionality  factor  of  1/R^  in  A 
is  the  integration  variable  whose  limits  are  both  on  the  same  side  of  zero. 
The  result  is 


Z2  " b y2  " 8 

- 2 log  (Uj  + y\  + x2)1/2  + y2)  (C-107) 

z2  " 0 y2  " 0 

which  is  unbounded  as  x^  approaches  zero,  that  is,  A ■*  21og  | x^ | ■* 

The  theoretical  consequences  of  the  above  discussion  are  very  interesting. 
Under  the  assumption  that  both  Jg  and  are  everywhere  integrable  in  either 
the  "t"  or  "a"  direction,  equation  (C-lOla)  shows  that  J may  not  be  Infinite 
except  at  a discontinuity  and  may  be  infinite  or  zero  but  not  finite  near 
a discontinuity.  The  latter  point  is  true  since  if  J(  la  assumed  to  be 
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bounded  away  from  zero  buc  not  Infinite  the  left-hand  side  is  finite  yet 
the  right-hand  side  displays  a logarithmic  singularity  as  it  approaches 
the  border.  The  former  point  follows  from  the  boundedness  of  the  integrals. 


To  continue  our  discussion  we  must  now  distinguish  between  junctions 
and  edges.  Our  discussion  of  the  behavior  of  a kernel  assumed  a junction 
and  showed  that  A was  negative.  For  an  edge  we  must  reverse  the  direction 
of  n causing  € to  reverse  direction  and  we  therefore  find  that  A is  posi- 
tive. This  distinction  is  critical,  for  by  rewriting  equation  (C-lOla) 


i Vr> 


K..  J (r  ) dSo  + 

A S O O i 


finite  terms 


(C-108) 


we  see  from  the  requirement  that  the  limit  of  Ja(r)  as  £ approaches  the 


border  be  the  same  if  we  approach  from  either  side  of  the  discontinuity 


that  when  A is  negative  Jg(rQ)  cannot  be  infinite,  but,  when  A is  positive 


J_(r_)  can  be  infinite.  Thus  J_  at  the  lunction  must  be  zero.  Furthermore, 
so  s 


the  need  for  the  singular  behavior  on  the  two  sides  of  equation  (C-108) 
to  exactly  balance,  greatly  restricts  the  type  of  singularity  that  can 
exist  at  the  edge.  This  restriction  enables  us  to  determine  the  nature  of 
the  singularity  as  will  be  illustrated  by  returning  to  the  wedge  problem 
considered  earlier.  We  assume  that  near  the  edge 


■ CI*2.IP  » Jg(*2^  ■ cz2’  -1  < P < o 


(C-109) 


(We  want  -1  < p so  that  J is  integrable  and  p < 0 so  that  J is  singular 

s s 


at  ■ 0.)  We  require 


„8  P+1  A A 

4"  J . J0  („;  ♦ ?2  ♦ ,2,3/2  2 ' l' 


(C-110) 


As  x^  -*■  0 then  x^/s  -*■  0,  x^/b  ■*  0 and  substitutions  y ■ yj/lxjJ  and 
z ■ Zj/lxjl  transform  the  integral 

00 

f zP+1  dz 

J 1 + z2 

o 

(C-lll) 

Since  -1  < p < 0 we  can  perform  contour  integration  for  the  final  integral 
obtaining 


c x. 


4ir 


00  ( 

// 


" zp+1  dz  dv  -l^ll1 
(1  + y2  + z2)3/2  2ir 


lx  lP 

t _ r 1 l1  3in[ir/2(p+l)  1 
A c 2 sin[ir(p  + 1)] 


(C-112) 


The  only  simultaneous  solution  to  equations  (C-112)  and  (C-110)  for  which 
-1  < p < 0 is  p - -1/3.  Thus  near  a right  angled  edge 

JS<*1>  * (C-113) 

This  result  may  be  derived  by  a number  of  different  approaches  (see 
for  example,  ref.  11).  Notice  that  when  c - 0 equation  (C-110)  is  auto- 
matically satisfied.  Whether  is  zero  or  infinite  for  an  edge  depends 
on  the  polarization  of  the  Incident  wave. 

A similar  discussion  for  equation  (C-lOlb)  shows  that  may  not 
become  unbounded  at  an  edge  since  for  an  edge  D is  negative.  To  eliminate 
the  possibility  of  an  unbounded  at  a Junction,  we  appeal  to  the  con- 
tinuity equation 


11.  Collin,  R.  E.,  Field  Theory  of  Guided  Waves.  McGraw-Hill  Book  Company, 
Inc.,  New  York,  1969,  pp.  18-20. 
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V • J_  “ iuxJ, 


(C-114) 


3J 
a 

3s 


iux? 


and  note  that  J , 3J  /3s,  and  o must  be  zero  for  a junction.  This  requires 
s s 

that  3Jt/3t  • 0 which  is  impossible  if  Jt(x^)  is  to  be  finite  for  |x^|  < 0 
and  infinite  at  » 0. 

The  numerical  consequences  of  the  above  discussions  are  very  straight- 
forward. A basic  premise  of  the  patch  zoning  technique  is  that  J can  be 

-1/3 

considered  to  be  nearly  constant  over  a patch.  If  J varies  as  z the 
, s 

patches  must  shrink  to  zero  area  in  order  to  keep  J nearly  constant  over 

s 

a patch;  this  would  require  an  infinite  number  of  zones  which  would  make 

costs  prohibitive  and  probably  make  the  matrix  very  ill-conditioned.  We 

must  therefore  accept  the  fact  that  calculated  values  for  J very  close 

s 

to  the  edge  will  be  relatively  Inaccurate.  Near  a junction  the  boundary 
condition  that  Jg  along  the  Intersection  be  zero  may  make  Jg  change 
rapidly  in  the  vicinity  of  the  intersection,  but  since  the  change  is  finite 
it  will  be  sufficient  to  slightly  Increase  the  zone  density  on  both  sides 
of  the  junction  to  obtain  substantially  accurate  results. 

The  presence  of  junctions  and  edges  has  profound  effects  on  our 
interpolation  procedure.  Since  an  interpolation  procedure  cannot  be  more 
reliable  than  the  data  it  has  to  work  with,  the  interpolation  will  fail 
near  an  endcap.  In  addition,  since  A has  a logarithmic  singularity  as  the 
observation  point  approaches  a discontinuity,  the  interpolation 
procedure  predicts  an  infinity  whether  the  discontinuity  is  an  edge  or 
a junction.  To  avoid  this  troublesome  fact  near  a junction  we  suggest 
that  no  attempt  be  made  to  determine  the  current  density  at  any  point 
closer  to  the  junction  than  the  reference  points  of  the  layer  zones 
closest  to  the  junction.  If  values  close  to  the  junction  are  desired 
the  reference  points  on  both  sides  of  the  junction  should  be  equally 
close  to  the  junction.  Also  a zone  should  not  be  significantly  smaller 
than  its  nearest  neighbors. 

In  concluding  this  section  we  would  like  to  dispell  any  doubts  that 
the  prediction  of  currents  away  from  the  endcaps  will  greatly  suffer  as 
a result  of  the  limited  accuracy  of  Jg  at  the  endcap.  The  Integral 
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equation  that  we  deal  with  is  very  stable.  The  kernel  is  singular  at 

r^  ■ r making  the  currents  in  the  innediate  vicinity  of  a point,  the 

dominant  Influence.  Therefore  the  singular  behavior  will  tend  to  damp  out 

as  we  move  away  from  the  endcap.  The  error  committed  by  approximating 

the  singularity  by  a finite  quantity  will  barely  effect  the  results.  As 

evidence  of  what  we  say,  we  presented  in  Section  IV,  the  results  of  an 

experimental  comparison  between  uncapped  cylinders,  flat  capped  cylinders 

and  hemispherically  capped  cylinders.  Even  though  the  theoretical  behavior 

-1/2  “1/3  0 

of  Jg  near  the  edge  is  x , x and  x for  the  three  respective  cases 
the  plots  of  Jt  do  not  show  any  appreciable  difference  due  to  the  varying 
boundary  conditions.  As  evidence  of  the  accuracy  of  our  programm  a quick 
look  at  figures  6 through  8 will  show  that  the  effects  of  the  external 
edge  did  not  greatly  effect  our  ability  to  predict  the  currents  on  a cylinder 
at  distances  as  close  to  the  edge  as  our  zone  size. 
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APPENDIX  D 


INTERSECTIONS 

In  this  appendix  we  consider  the  intersection  between  the  elliptical 
bodies-components  of  the  aircraft  and  show  how  their  presence  modifies 
(a)  the  coordinates  at  the  centers  of  zones  adjacent  to  the  intersections 

and  (b)  the  matrix  element  associated  with  these  zones. 

1.  COORDINATES 

a.  Zones  Adjacent  to  Fuselage  - Wing  Intersections 
(1)  Zones  on  Fuselage  (figure  Dl) 

We  assume  that  the  planes  x » x^  and  z **  0 form  walls  for  zones  on 

the  fuselage.  Also  only  one  layer  of  zones  is  intersected  for  both  z > 0 

and  z < 0.  (That  is,  zones  that  would  exist  in  the  absence  of  the  inter- 
section.) If  the  center  of  a zone  (as  defined  in  Appendix  C,  equation 
(C-18),  where  no  intersections  were  assumed)  has  an  x-coordinate  such  that 
x < x02  “ a2  or  x > XQ2  + a2  tben  coordinates  at  the  center  of  this 
zone  are  unaffected  by  the  presence  of  the  Intersection. 


x02  + a2  > x > x02  " 

a 2 then: 

x(a)  - Xj^ 

8x(°)  " 0 

t (a) 

X 

b,  cos$ 

y(ct)  ■ b^  sin$ 

8y(a)  “ nxW) 

ty(a) 

a sin<f> 

z(a)  ■ -a^  cos<|» 

•«<0)  ’ «.(♦) 

tz(a) 

■ i ~2  <P  (+  for  z > 0 and  - for  z < 0) 


* ' 810-1  i"  z*  “ b2(1  " *T 
1 \ *2 


X - x„„  - 


02  " *k 
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WING 


and  <J>j , x^,  N^(<{>)  are  given  by  equations  (C-18) . 

(2)  Zones  on  Wing  (figure  Dl) 

The  new  coordinates  at  the  center  of  the  zones  adjacent  to  the 
intersection  are 


x(ot)  - Xq2  + a2  cos<^ 


/■\  1 * 

y(°0  - yk  - 2 y 


z(a)  ■ -b2  sin$j 


-a  sin4>< 

• Cx<a>  - 0 


sy(a)  » 0 


, ty(a)  - 1 


b»  cos<f) 

3 _(a)  „ ^ , t (a)  - 0 

z (<Pj  ) z 


where  y^»  N2 (♦)  are  given  by  equations  (C-22)  and 


* [ / b2  Sin2<* 
1 ■‘•r  "t 


n 


b.  Zones  Adjacent  to  Fuselage  - Horizontal  Stabilizer 

(1)  Zones  on  Fuselage 

The  new  coordinates  of  the  zones  are  derived  by  the  same  assumptions 
as  for  the  wing  and  are  identical  to  the  ones  given  by  equations  (D-l) 
provided  we  replace  Xq2>  b2>  a2  by  x^,  b^,  a^  respectively. 

(2)  Zones  on  Horizontal  Stabilizer 

The  new  coordinates  are  derived  by  the  same  assumptions  as  for  the 
wing  and  are  given  by  equations  (D-2)  provided  we  replace  xQ2>  a2>  b2, 

N2  by  Xq^  bj.  respectively  and  refer  to  equations  (C-24)  instead 

of  (C-22). 

c.  Zones  Adjacent  to  Fuselage  - Vertical  Stabilizer  Intersection 

(1)  Zones  on  Fuselage  (figure  D2) 


We  assume  that  the  y ■ 0 and  x - x^  planes  form  natural  boundaries 
for  zones  on  the  fuselage.  As  with  the  fuselage-wing  intersection  we 
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VERTICAL  STABILIZER 


again  assume  there  is  only  one  layer  of  zones  intersected.  The  new 
coordinates  are  given  by 


x(a) 

' \ 

3 (a) 

X 

- o 

tx(a) 

cos$ 

y(a) 

■ a^  sin<|) 

sy(o) 

\<A)  * 

ty(a) 

b^  sin<(> 

z(a) 

■ -b^  cos<f> 

s2(°t) 

N.(o)  ’ 

tz(a) 

where 


1 * 
' i* 


♦ - sin-1  *T  y*  - b4(l  - 


a/2 


X * X04  - ^ 

and  xk’  are  given  by  equations  (C-18). 

(2)  Zones  on  Vertical  Stabilizer  (figure  D2) 
The  new  coordinates  are 


i I 

U 5 


r * 


x(0t)  " x04  + a4  C08<^j 


y(a)  - b^  sin^ 


a sin0 

8x(a) 5^)  • Cx<a)  ■ 1 

b,  cos<}> 

Sy(0t)  “ N4(^)J  ’ ty(a)  " C 


1 * 


Z(a)  “ 2k  ‘ 2 Z 


s (a)  - 0 
z 


, tz(o)  - ] 


where 


* 

z - a. 


1 - 1 - 


2 2 
bf  sin \4 
A Li 


a/2n 


end  ♦j • z^,  are  given  by  equation  (C-28) , 
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(D-3) 


(D-4) 


2.  MATRIX  ELEMENTS 


As  we  explained  In  Appendix  C we  calculated  the  matrix  elements 
(equations  (C-63)  through  (C-96))  by  employing  auxiliary  local  coordinate 
systems,  attached  to  each  cylindrical  body,  such  that  x^  was  along  the 
axis  of  the  body.  The  formulas  giving  the  matrix  elements  must  eventually 
be  expressed  in  terms  of  the  global  coordinate  system  (figure  1)  through 
transformation  (C-30a)  and  (C-30b).  (In  the  computer  code  the  matrix 
elements  are  calculated  in  the  global  system.)  In  this  subsection,  we 
present  the  matrix  elements  in  terms  of  the  local  coordinate  systems. 

Notice  that  only  zones  on  the  walls  are  affected. 

a.  Zones  Adjacent  to  Fuselage  - Wing  Intersection 
(1)  Zones  on  Fuselage  (figure  D3) 

Figure  21  shows  how  we  subdivide  the  zones  adjacent  to  the  intersection 
into  subzones  (1)  and  (2).  First  we  will  restrict  our  attention  to  z > 0, 
x < Xq2  (figure  21). 

(a)  Region  z > 0,  x < xQ2 

The  matrix  elements  A*  are  given  by  equations  (C-63)  through  (C-65) 
with  the  following  changes.  For  subzones  labeled  (1)  in  figure  D3  we 
replace  <j>^,  $2  (as  they  appear  in  these  equations)  by  4> , $2  respectively 
where 


and  $2  defines  a straight  boundary  of  the  zone.  For  subzones  labeled 
(2)  in  figure  21,  <p2,  x1Q1,  x1Q2  (as  they  appear  in  equations  (C-63) 

through  (C-65))  are  replaced  by  <pj,  x10l’  (figure  D3)  respec- 

tively, where  for  xQ2  - x1Q1  <_  a2. 
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-> ... 


4 

tcj 

3 


►i  - cos'1  if  f1 


(x02  " ’‘lOl* 


(2  2 
a.  cos  4> 

1 - - b2 


(D-6) 


Notice  that  when  > a^,  and  define  the  two  straight 

^-boundaries  of  the  leftmost  zone  in  figure  D3. 

(b)  Region  z < 0,  x < Xq2  (figure  D4) 

Keeping  in  mind  that  in  q < 0 we  still  have  < <p7  < <p~,  the  matrix 
element  A"  are  given  by  equations  (C-63)  through  (C-65)  with  the  following 
changes.  For  subzones  (1)  in  figure  D4  we  replace  <J>^,  (as  they  appear 
in  these  equations)  by  where 


-i°'1  i 1 


(x02  ~ X102) 


(D-7) 


and  defines  a straight  boundary.  For  subzones  labeled  (2)  in  figure  D4 

$1*  $2*  xioi*  xl02  ^aa  they  aPPear  in  equations  (C-63)  through  (C-65) 
are  replaced  by  0^.  <J>2*  xioi*  (figure  D4)  respectively,  where  for 


x02  ~ x101  - a2 


•u-1 

. 1 V *2 


X1  (^)  “ X02  “ ®2  ( 1 " 


2 2 
a cos  $ 


(D-8) 


I 
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When  Xq2  ~ x^gi  > a2*  ^1  an<*  ^2  ^e^ine  c^e  two  straight  ^-boundaries 
of  the  leftmost  zone  in  figure  D4  for  z < 0. 

(c)  Region  x > xQ2,  z > 0 

The  matrix  elements  are  derived  by  a procedure  identical  to  that 
employed  in  region  x < x^,  z > 0 provided  we  interchange  the  definitions 
for  <J>^  and  *3  (that  is,  *3  is  now  defined  by  equation  (D-5)  and  *3  by 
(D-6))  and  define 


x1(4>) 


(D-9) 


(d)  Region  x > xQ2,  z < 0 

The  results  are  identical  to  those  in  region  x < Xq2>  z < 0 provided 
we  interchange  the  definitions  of  *2  and  *3  (that  is,  *2  is  now  defined 
by  equation  (D-7)  and  *3  by  (D-8))  and  define  x^(*)  by  (D-9). 

+ ■*•  + + 

Matrix  elements  B~,  C~,  and  D“  can  be  calculated  similarly  as  A" 

provided  we  use  the  corresponding  equations  given  in  Appendix  C. 


(2)  Zones  on  Wing  (figure  D3) 


Matrix  elements  A are  given  by  equations  (C-71)  through  (C-73) 
provided  we  replace  *102  (in  these  equations)  by  x^*),  x1Q2 
(figure  D3)  respectively,  where 


XjOfr)  » a 


1 


b2  sin2*  \1/2 
) 


(D-10) 


Matrix  elements  B~ , C~  and  D~  can  be  obtained  similarly  provided  we 
use  the  corresponding  equations  given  in  Appendix  C. 

b.  Zones  Adjacent  to  Fuselage  - Horizontal  Stabilizer  Interaction 
Matrix  elements  A1,  B± , C~  and  D*  for  zones  on  both  the  fuselage 
and  the  horizontal  stabilizer  are  calculated  in  a manner  identical  to 
the  one  employed  for  zones  adjacent  to  the  fuaelage-wlng  interaction. 
Naturally,  we  must  replace  xQ2,  a2,  b2  by  xQ3,  a3,  b3  respectively. 
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c.  Zones  Adjacent  to  Fuselage  - Vertical  Stabilizer  Intersection 

(1)  Zones  on  Fuselage  (figure  D5) 

The  matrix  element  A-  are  given  by  equations  (C— 63)  through  (C-65) 
with  the  following  changes.  For  subzones  labeled  (1)  in  figure  05  we 
replace  <P2  (as  they  appear  In  the  equations)  by  <t>^,  <P2  (figure  D5) 
respectively  where 


(0-11) 


and  $2  defines  the  other  boundary  of  the  zone.  For  subzones  labeled  (2) 
in  figure  05 , we  replace  <J>2  (as  they  appear  in  equations  (C-63) 
through  (C-65)  by  (figure  D5)  respectively  and  x1Q1>  x1Q2  (as 

they  appear  in  equations  (C-63)  through  (C-65))  by  x1Q1  (figure  D5), 
x^(<t>)  respectively  where 


2 2 
b*  sln^ 


1/2 


*!<$>  - xQ4  * a4|l -7j j , (-)  for  x < xQ4,  (+)  for  x > xQ4 


- sin 


-1 


(x04  ~ *101^ 


1/2! 


' x < x04’  x04  “ x101  < a4 


• sin 


-1 


^4  / (XQ4  ~ x102)< 

.1 V t 


,1/2! 


x > x04*  x04  “ x101  < a4 


(D-12) 


VERTICAL  STABILIZER 


Notice  that  when  xQ4  - x1Q1  > a4  in  region  x < xQ4  or  x1Q2  - xQ4  > a4 
in  region  x > x^4,  <J>^  and  <|>2  define  the  two  straight  (^-boundaries  for 
the  leftmost  zone  in  region  x < xQ4  or  the  rightmost  zone  in  region 
x > xQ4  (figure  D5) . 

+ + + 

Matrix  elements  B~,  C”,  and  D~  can  be  derived  similarly  by  using 
the  same  substitution  as  for  A-  in  the  relevant  equation  given  in 
Appendix  C. 

(2)  Zones  on  Vertical  Stabilizer  (figure  D5) 

Matrix  elements  A-  are  given  by  equations  (C-78)  through  (C-80) 
provided  we  replace  *101>  x1Q2  (in  these  equations)  by  x1(4>),  x1Q2 
(figure  D5)  respectively,  where 


b 2 sin2<( 

Xl(*}  " "l  l1 ^2 


1/2 


(D— 1 3 ) 


Matrix  elements  B~,  C~  and  D~  can  be  derived  similarly  as  we  remarked 
right  after  equation  (D-ll). 
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